kNN algorithm improved. IO Defined
authorabercuci <abercuci@f7af4fe6-9843-0410-8265-dc069ae4e863>
Fri, 21 Sep 2007 15:06:13 +0000 (15:06 +0000)
committerabercuci <abercuci@f7af4fe6-9843-0410-8265-dc069ae4e863>
Fri, 21 Sep 2007 15:06:13 +0000 (15:06 +0000)
STAT/TKDInterpolator.cxx
STAT/TKDInterpolator.h
STAT/TKDTree.cxx
STAT/TKDTree.h

index 9be1a2f..ad7f08d 100644 (file)
@@ -13,9 +13,8 @@
 #include "TRandom.h"
 #include "TROOT.h"
 
-
-
 ClassImp(TKDInterpolator)
+ClassImp(TKDInterpolator::TKDNodeInfo)
 
 /////////////////////////////////////////////////////////////////////
 // Memory setup of protected data memebers
@@ -24,6 +23,8 @@ ClassImp(TKDInterpolator)
 //
 // fRefValues : evaluation value/error of PDF for each terminal node of underlying KD Tree.
 // | 1st terminal node (value) | 2nd terminal node (value) | ... | 1st terminal node (error) | 2nd terminal node (error) | ...
+//
+// status = |0|0|0|0|0|1(tri-cubic weights)|1(STORE)|1 INT(0 COG )|
 /////////////////////////////////////////////////////////////////////
 
 //_________________________________________________________________
@@ -31,8 +32,13 @@ TKDInterpolator::TKDInterpolator() : TKDTreeIF()
        ,fNTNodes(0)
        ,fRefPoints(0x0)
        ,fRefValues(0x0)
+       ,fCov(0x0)
+       ,fPar(0x0)
+       ,fPDFstatus(0x0)
+       ,fStatus(4)
+       ,fLambda(0)
        ,fDepth(-1)
-       ,fTmpPoint(0x0)
+       ,fBuffer(0x0)
        ,fKDhelper(0x0)
        ,fFitter(0x0)
 {
@@ -45,8 +51,13 @@ TKDInterpolator::TKDInterpolator(Int_t npoints, Int_t ndim, UInt_t bsize, Float_
        ,fNTNodes(GetNTerminalNodes())
        ,fRefPoints(0x0)
        ,fRefValues(0x0)
+       ,fCov(0x0)
+       ,fPar(0x0)
+       ,fPDFstatus(0x0)
+       ,fStatus(4)
+       ,fLambda(0)
        ,fDepth(-1)
-       ,fTmpPoint(0x0)
+       ,fBuffer(0x0)
        ,fKDhelper(0x0)
        ,fFitter(0x0)
 {
@@ -57,12 +68,17 @@ TKDInterpolator::TKDInterpolator(Int_t npoints, Int_t ndim, UInt_t bsize, Float_
 
 
 //_________________________________________________________________
-TKDInterpolator::TKDInterpolator(TTree *t, const Char_t *var, const Char_t *cut, UInt_t bsize) : TKDTreeIF()
+TKDInterpolator::TKDInterpolator(TTree *t, const Char_t *var, const Char_t *cut, UInt_t bsize, Long64_t nentries, Long64_t firstentry) : TKDTreeIF()
        ,fNTNodes(0)
        ,fRefPoints(0x0)
        ,fRefValues(0x0)
+       ,fCov(0x0)
+       ,fPar(0x0)
+       ,fPDFstatus(0x0)
+       ,fStatus(4)
+       ,fLambda(0)
        ,fDepth(-1)
-       ,fTmpPoint(0x0)
+       ,fBuffer(0x0)
        ,fKDhelper(0x0)
        ,fFitter(0x0)
 {
@@ -73,23 +89,25 @@ TKDInterpolator::TKDInterpolator(TTree *t, const Char_t *var, const Char_t *cut,
 // 
 
        TObjArray *vars = TString(var).Tokenize(":");
-       fNDim = vars->GetEntriesFast();
+       fNDim = vars->GetEntriesFast(); fNDimm = 2*fNDim;
        if(fNDim > 6/*kDimMax*/) Warning("TKDInterpolator(TTree*, const Char_t, const Char_t, UInt_t)", Form("Variable number exceed maximum dimension %d. Results are unpredictable.", 6/*kDimMax*/));
        fBucketSize = bsize;
 
        Int_t np;
        Double_t *v;
        for(int idim=0; idim<fNDim; idim++){
-               if(!(np = t->Draw(((TObjString*)(*vars)[idim])->GetName(), cut, "goff"))){
-                       Warning("TKDInterpolator(TTree*, const Char_t, const Char_t, UInt_t)", Form("Can not access data for %s", ((TObjString*)(*vars)[idim])->GetName() ));
+               if(!(np = t->Draw(((TObjString*)(*vars)[idim])->GetName(), cut, "goff", nentries, firstentry))){
+                       Warning("TKDInterpolator(TTree*, const Char_t, const Char_t, UInt_t)", Form("Can not access data for keys %s. Key defined on tree :", ((TObjString*)(*vars)[idim])->GetName() ));
+                       TIterator *it = (t->GetListOfLeaves())->MakeIterator();
+                       TObject *o;
+                       while(o = (*it)()) printf("\t%s\n", o->GetName());
                        continue;
                }
                if(!fNpoints){
                        fNpoints = np;
                        Info("TKDInterpolator(TTree*, const Char_t, const Char_t, UInt_t)", Form("Allocating %d data points in %d dimensions.", fNpoints, fNDim));
-                       //Float_t *mem = new Float_t[fNDim*fNpoints];
                        fData = new Float_t*[fNDim];
-                       for(int idim=0; idim<fNDim; idim++) fData[idim] = new Float_t[fNpoints]; //&mem[idim*fNpoints];
+                       for(int idim=0; idim<fNDim; idim++) fData[idim] = new Float_t[fNpoints];
                        kDataOwner = kTRUE;
                }
                v = t->GetV1();
@@ -103,9 +121,15 @@ TKDInterpolator::TKDInterpolator(TTree *t, const Char_t *var, const Char_t *cut,
 //_________________________________________________________________
 TKDInterpolator::~TKDInterpolator()
 {
+       if(fCov){
+               delete [] fPar;
+               delete [] fCov;
+               delete [] fPDFstatus;
+       }
+       
        if(fFitter) delete fFitter;
        if(fKDhelper) delete fKDhelper;
-       if(fTmpPoint) delete [] fTmpPoint;
+       if(fBuffer) delete [] fBuffer;
        
        if(fRefPoints){
                for(int idim=0; idim<fNDim; idim++) delete [] fRefPoints[idim] ;
@@ -122,7 +146,8 @@ void TKDInterpolator::Build()
 //  - corresponding PDF values
 
        if(!fBoundaries) MakeBoundaries();
-       
+       fLambda = 1 + fNDim + fNDim*(fNDim+1)/2;
+
        // allocate memory for data
        fRefValues = new Float_t[fNTNodes];
        fRefPoints = new Float_t*[fNDim];
@@ -162,76 +187,217 @@ void TKDInterpolator::Build()
        }
 }
 
+//__________________________________________________________________
+void TKDInterpolator::GetStatus()
+{
+       printf("Interpolator Status :\n");
+       printf("  Method : %s\n", fStatus&1 ? "INT" : "COG");
+       printf("  Store  : %s\n", fStatus&2 ? "YES" : "NO");
+       printf("  Weights: %s\n", fStatus&4 ? "YES" : "NO");
+
+       printf("nodes %d\n", fNTNodes);        //Number of evaluation data points
+       printf("refs 0x%x\n", fRefPoints);    //[fNDim][fNTNodes]
+       printf("vals 0x%x\n", fRefValues);     //[fNTNodes]
+       for(int i=0; i<fNTNodes; i++){
+               for(int idim=0; idim<fNDim; idim++) printf("%f ", fRefPoints[idim][i]);
+               printf("[%f]\n", fRefValues[i]);
+       }
+
+       printf("cov 0x%x\n", fCov);           //[fNTNodes] cov matrix array for nodes
+       printf("pars 0x%x\n", fPar);           //[fNTNodes] parameters array for nodes
+       for(int i=0; i<fNTNodes; i++){
+               printf("%d ", i);
+               for(int ip=0; ip<3; ip++) printf("p%d[%f] ", ip, fPar[i](ip));
+               printf("\n");
+       }
+       printf("pdf 0x%x\n", fPDFstatus);     //[fNTNodes] status bit for node's PDF
+}
+
 //_________________________________________________________________
-Double_t TKDInterpolator::Eval(const Double_t *point, Int_t npoints)
+Double_t TKDInterpolator::Eval(const Double_t *point, Double_t &result, Double_t &error)
 {
-// Evaluate PDF at k-dimensional position "point". The initial number of
-// neighbour estimation points is set to "npoints"
+// Evaluate PDF for "point". The result is returned in "result" and error in "error". The function returns the chi2 of the fit.
+//
+// Observations:
+//
+// 1. The default method used for interpolation is kCOG.
+// 2. The initial number of neighbors used for the estimation is set to Int(alpha*fLambda) (alpha = 1.5)
+                       
+       Float_t pointF[50]; // local Float_t conversion for "point"
+       for(int idim=0; idim<fNDim; idim++) pointF[idim] = (Float_t)point[idim];
+       Int_t node = FindNode(pointF) - fNnodes;
+       if(fPDFstatus && (fStatus&1) && fPDFstatus[node]) return CookPDF(point, node, result, error);
+
+       // Allocate memory
+       if(!fBuffer) fBuffer = new Double_t[2*fLambda];
+       if(!fKDhelper) fKDhelper = new TKDTreeIF(fNTNodes, fNDim, 30, fRefPoints);
+       if(!fFitter) SetIntInterpolation(kFALSE);
        
+       // generate parabolic for nD
+       //Float_t alpha = Float_t(2*lambda + 1) / fNTNodes; // the bandwidth or smoothing parameter
        //Int_t npoints = Int_t(alpha * fNTNodes);
        //printf("Params : %d NPoints %d\n", lambda, npoints);
        // prepare workers
-       if(!fTmpPoint) fTmpPoint = new Double_t[fNDim];
-       if(!fKDhelper) fKDhelper = new TKDTreeIF(GetNTerminalNodes(), fNDim, npoints, fRefPoints);
-       if(!fFitter){
-               // generate parabolic for nD
-               
-               // calculate number of parameters in the parabolic expresion
-               Int_t lambda = 1 + fNDim + fNDim*(fNDim+1)/2;
-               //Float_t alpha = Float_t(2*lambda + 1) / fNTNodes; // the bandwidth or smoothing parameter
-               TString formula("1");
-               for(int idim=0; idim<fNDim; idim++){
-                       formula += Form("++x[%d]", idim);
-                       for(int jdim=idim; jdim<fNDim; jdim++) formula += Form("++x[%d]*x[%d]", idim, jdim);
-               }
-               fFitter = new TLinearFitter(lambda, formula.Data());
-               Info("Eval(const Double_t*, Int_t)", Form("Using %s for local interpolation.", formula.Data()));
-       }
 
-       Float_t pointF[50];
-       for(int idim=0; idim<fNDim; idim++) pointF[idim] = point[idim];
-       Int_t istart = 0;
-       Int_t *index;
-       Float_t dist, d0, w0, w;
-       Double_t uncertainty = TMath::Sqrt(1./fBucketSize); 
-       fFitter->ClearPoints();
+       Int_t *index,  // indexes of NN 
+             ipar,    // local looping variable
+                               npoints = Int_t(1.5*fLambda); // number of data points used for interpolation
+       Float_t *dist, // distances of NN
+                                       d,     // NN normalized distance
+                                       w0,    // work
+                                       w;     // tri-cubic weight function
+       Double_t sig   // bucket error 
+               = TMath::Sqrt(1./fBucketSize);
        do{
+               // find nearest neighbors
+               for(int idim=0; idim<fNDim; idim++) pointF[idim] = (Float_t)point[idim];
                if(!fKDhelper->FindNearestNeighbors(pointF, npoints+1, index, dist)){
                        Error("Eval()", Form("Failed retriving %d neighbours for point:", npoints));
                        for(int idim=0; idim<fNDim; idim++) printf("%f ", point[idim]);
                        printf("\n");
                        return -1;
                }
-               for(int in=istart; in<npoints; in++){
-                       //printf("%d index[%2d] x(", in, index[in]);
-                       d0 = 0.;
-                       for(int idim=0; idim<fNDim; idim++){
-                               fTmpPoint[idim] = fRefPoints[idim][index[in]];
-                               //printf("%6.4f ", fTmpPoint[idim]);
-                               d0 += TMath::Abs(fTmpPoint[idim] - point[idim]);
+               // add points to fitter
+               fFitter->ClearPoints();
+               for(int in=0; in<npoints; in++){
+                       if(fStatus&1){ // INT
+                               for(int idim=0; idim<fNDim; idim++) pointF[idim] = fRefPoints[idim][index[in]];
+                               Float_t *bounds = GetBoundary(FindNode(pointF));
+                               
+                               ipar = 0;
+                               for(int idim=0; idim<fNDim; idim++){
+                                       fBuffer[ipar++] = .5*(bounds[2*idim] + bounds[2*idim+1]);
+                                       fBuffer[ipar++] = (bounds[2*idim]*bounds[2*idim] + bounds[2*idim] * bounds[2*idim+1] + bounds[2*idim+1] * bounds[2*idim+1])/3.;
+                                       for(int jdim=idim+1; jdim<fNDim; jdim++) fBuffer[ipar++] = (bounds[2*idim] + bounds[2*idim+1]) * (bounds[2*jdim] + bounds[2*jdim+1]) * .25;
+                               }
+                       } else { // COG
+                               for(int idim=0; idim<fNDim; idim++) fBuffer[idim] = fRefPoints[idim][index[in]];
                        }
-                       d0 /= dist;
-                       w0 = (1. - d0*d0*d0);
-                       w = w0*w0*w0;
 
-                       //printf(") f = %f [%f] d0 = %6.4f   w = %6.4f  \n", fRefValues[index[in]], TMath::Log(fRefValues[index[in]]), d0, w);
-                       fFitter->AddPoint(fTmpPoint, TMath::Log(fRefValues[index[in]]), uncertainty/w);
+                       // calculate tri-cubic weighting function
+                       if(fStatus&4){
+                               d = dist[in]/ dist[npoints];
+                               w0 = (1. - d*d*d); w = w0*w0*w0;
+                       } else w = 1.;
+                        
+                       //for(int idim=0; idim<fNDim; idim++) printf("%f ", fBuffer[idim]);
+                       //printf("\nd[%f] w[%f] sig[%f]\n", d, w, sig);
+                       fFitter->AddPoint(fBuffer, fRefValues[index[in]], fRefValues[index[in]]*sig/w);
                }
-               istart = npoints;
                npoints += 4;
        } while(fFitter->Eval());
 
-       // calculate evaluation
-       Int_t ipar = 0;
-       Double_t result = fFitter->GetParameter(ipar++);
+       // retrive fitter results
+       TMatrixD cov(fLambda, fLambda);
+       TVectorD par(fLambda);
+       fFitter->GetCovarianceMatrix(cov);
+       fFitter->GetParameters(par);
+       Double_t chi2 = fFitter->GetChisquare()/(npoints - 4 - fLambda);
+
+       // store results
+       if(fStatus&2 && fStatus&1){
+               fPar[node] = par;
+               fCov[node] = cov;
+               fPDFstatus[node] = 1;
+       }
+               
+       // Build df/dpi|x values
+       Double_t *fdfdp = &fBuffer[fLambda];
+       ipar = 0;
+       fdfdp[ipar++] = 1.;
        for(int idim=0; idim<fNDim; idim++){
-               result += fFitter->GetParameter(ipar++)*point[idim];
-               for(int jdim=idim; jdim<fNDim; jdim++) result += fFitter->GetParameter(ipar++)*point[idim]*point[jdim];
+               fdfdp[ipar++] = point[idim];
+               for(int jdim=idim; jdim<fNDim; jdim++) fdfdp[ipar++] = point[idim]*point[jdim];
        }
-       //printf("\tResult : %f [%f]\n", TMath::Exp(result), result);
-       return TMath::Exp(result);
+
+       // calculate estimation
+       result =0.; error = 0.;
+       for(int i=0; i<fLambda; i++){
+               result += fdfdp[i]*par(i);
+               for(int j=0; j<fLambda; j++) error += fdfdp[i]*fdfdp[j]*cov(i,j);
+       }       
+       error = TMath::Sqrt(error);
+
+       return chi2;
 }
 
+// //_________________________________________________________________
+// Double_t TKDInterpolator::Eval1(const Double_t *point, Int_t npoints, Double_t &result, Double_t &error)
+// {
+// // Evaluate PDF at k-dimensional position "point". The initial number of
+// // neighbour estimation points is set to "npoints". The default method
+// // used for interpolation is kCOG.
+// 
+//     // calculate number of parameters in the parabolic expresion
+//     Int_t lambda = 1 + fNDim + fNDim*(fNDim+1)/2;
+// 
+//     if(!fBuffer) fBuffer = new Double_t[lambda-1];
+//     if(!fKDhelper) fKDhelper = new TKDTreeIF(GetNTerminalNodes(), fNDim, npoints, fRefPoints);
+// 
+//     if(!fFitter) fFitter = new TLinearFitter(lambda, Form("hyp%d", fNDim+1));
+//     else fFitter->SetFormula(Form("hyp%d", fNDim+1));
+// 
+// 
+//     Float_t pointF[50];
+//     for(int idim=0; idim<fNDim; idim++) pointF[idim] = point[idim];
+//     Int_t istart = 0;
+//     Int_t *index, ipar;
+//     Float_t *bounds, *dist, *w = new Float_t[fNDim];
+//     Double_t uncertainty = TMath::Sqrt(1./fBucketSize);
+//     fFitter->ClearPoints();
+//     do{
+//             if(!fKDhelper->FindNearestNeighbors(pointF, npoints+1, index, dist)){
+//                     Error("Eval()", Form("Failed retriving %d neighbours for point:", npoints));
+//                     for(int idim=0; idim<fNDim; idim++) printf("%f ", point[idim]);
+//                     printf("\n");
+//                     return -1;
+//             }
+//             for(int in=istart; in<npoints; in++){
+//                     for(int idim=0; idim<fNDim; idim++) w[idim] = fRefPoints[idim][index[in]];
+//                     bounds = GetBoundary(FindNode(w));
+// 
+//                     ipar = 0;
+//                     for(int idim=0; idim<fNDim; idim++){
+//                             fBuffer[ipar++] = .5*(bounds[2*idim] + bounds[2*idim+1]);
+//                             fBuffer[ipar++] = (bounds[2*idim]*bounds[2*idim] + bounds[2*idim] * bounds[2*idim+1] + bounds[2*idim+1] * bounds[2*idim+1])/3.;
+//                             for(int jdim=idim+1; jdim<fNDim; jdim++) fBuffer[ipar++] = (bounds[2*idim] + bounds[2*idim+1]) * (bounds[2*jdim] + bounds[2*jdim+1]) * .25;
+//                     }
+// 
+//                     fFitter->AddPoint(fBuffer, fRefValues[index[in]], fRefValues[index[in]]*uncertainty);
+//             }
+//             istart = npoints;
+//             npoints += 4;
+//     } while(fFitter->Eval());
+//     delete [] w;
+// 
+//     // calculate evaluation
+// //  fFitter->PrintResults(3);
+//     TMatrixD cov(lambda, lambda);
+//     TVectorD par(lambda);
+//     fFitter->GetCovarianceMatrix(cov);
+//     fFitter->GetParameters(par);
+// 
+//     // Build temporary array to keep values df/dpi|x
+//     Double_t f[100];
+//     ipar = 0;
+//     f[ipar++] = 1.;
+//     for(int idim=0; idim<fNDim; idim++){
+//             f[ipar++] = point[idim];
+//             for(int jdim=idim; jdim<fNDim; jdim++) f[ipar++] = point[idim]*point[jdim];
+//     }
+//     result =0.; error = 0.;
+//     for(int i=0; i<lambda; i++){
+//             result += f[i]*par[i];
+//             for(int j=0; j<lambda; j++) error += f[i]*f[j]*cov(i,j);
+//     }
+//     error = TMath::Sqrt(error);
+//     Double_t chi2 = fFitter->GetChisquare()/(npoints - 4 - lambda);
+// 
+//     for(int ipar=0; ipar<lambda; ipar++) printf("%d %8.6e %8.6e\n", ipar, par[ipar], TMath::Sqrt(cov(ipar, ipar)));
+//     printf("result %6.3f +- %6.3f [%f]\n", result, error, chi2);
+//     return chi2;
+// }
+
 
 //_________________________________________________________________
 void TKDInterpolator::DrawNodes(UInt_t ax1, UInt_t ax2, Int_t depth)
@@ -335,3 +501,88 @@ void TKDInterpolator::DrawNode(Int_t tnode, UInt_t ax1, UInt_t ax2)
        return;
 }
 
+
+//__________________________________________________________________
+void TKDInterpolator::SetIntInterpolation(const Bool_t on)
+{
+// Set interpolation bit to "on" and build/delete memory
+       
+       if(on) fStatus += fStatus&1 ? 0 : 1;
+       else fStatus += fStatus&1 ? -1 : 0;
+       TString formula;
+       if(on) formula = Form("hyp%d", fLambda-1);
+       else {
+               formula = "1";
+               for(int idim=0; idim<fNDim; idim++){
+                       formula += Form("++x[%d]", idim);
+                       for(int jdim=idim; jdim<fNDim; jdim++) formula += Form("++x[%d]*x[%d]", idim, jdim);
+               }
+       }
+       if(!fFitter) fFitter = new TLinearFitter(fLambda, formula.Data());
+       else fFitter->SetFormula(formula.Data());
+}
+
+
+//_________________________________________________________________
+void TKDInterpolator::SetSetStore(const Bool_t on)
+{
+// Set store bit to "on" and build/delete memory
+       
+       if(on){
+               fStatus += fStatus&2 ? 0 : 2;
+               if(!fCov){
+                       fPDFstatus = new Bool_t[fNTNodes];
+                       fCov = new TMatrixD[fNTNodes];
+                       fPar = new TVectorD[fNTNodes];
+                       for(int i=0; i<fNTNodes; i++){
+                               fPDFstatus[i] = kFALSE;
+                               fCov[i].ResizeTo(fLambda, fLambda);
+                               fPar[i].ResizeTo(fLambda);
+                       }
+               }
+       } else {
+               fStatus += fStatus&2 ? -2 : 0;
+               if(fCov){
+                       delete [] fPar;
+                       delete [] fCov;
+                       delete [] fPDFstatus;
+               }
+       }
+}
+
+//_________________________________________________________________
+void TKDInterpolator::SetUseWeights(const Bool_t on)
+{
+       if(on) fStatus += fStatus&4 ? 0 : 4;
+       else fStatus += fStatus&4 ? -4 : 0;
+}
+
+
+//_________________________________________________________________
+Double_t TKDInterpolator::CookPDF(const Double_t *point, const Int_t node, Double_t &result, Double_t &error)
+{
+// Recalculate the PDF for one node from the results of interpolation (parameters and covariance matrix)
+
+       Info("CookPDF()", Form("Called for node %d", node));
+
+       if(!fBuffer) fBuffer = new Double_t[2*fLambda];
+       Double_t *fdfdp = &fBuffer[fLambda];
+       Int_t ipar = 0;
+       fdfdp[ipar++] = 1.;
+       for(int idim=0; idim<fNDim; idim++){
+               fdfdp[ipar++] = point[idim];
+               for(int jdim=idim; jdim<fNDim; jdim++) fdfdp[ipar++] = point[idim]*point[jdim];
+       }
+
+       // calculate estimation
+       result =0.; error = 0.;
+       for(int i=0; i<fLambda; i++){
+               result += fdfdp[i]*fPar[node](i);
+               for(int j=0; j<fLambda; j++) error += fdfdp[i]*fdfdp[j]*fCov[node](i,j);
+       }       
+       error = TMath::Sqrt(error);
+       printf("result[CookPDF] %6.3f +- %6.3f\n", result, error);
+
+       return 0.;
+}
+
index e1d0ab3..145b20c 100644 (file)
@@ -4,39 +4,66 @@
 #ifndef ROOT_TKDTree
 #include "TKDTree.h"
 #endif
+#ifndef ROOT_TVectorD
+#include "TVectorD.h"
+#endif
+#ifndef ROOT_TMatrixD
+#include "TMatrixD.h"
+#endif
+
 
 class TTree;
 class TLinearFitter;
 class TKDInterpolator : public TKDTreeIF
 {
 public:
+       struct TKDNodeInfo
+       {
+
+               ClassDef(TKDNodeInfo, 1) // node info for interpolator
+       };
+public:
        TKDInterpolator();
-       TKDInterpolator(TTree *t, const Char_t *var, const Char_t *cut = 0, UInt_t bsize = 100);
+       TKDInterpolator(TTree *t, const Char_t *var, const Char_t *cut = 0, UInt_t bsize = 100, Long64_t nentries = 1000000000, Long64_t firstentry = 0);
        TKDInterpolator(Int_t npoints, Int_t ndim, UInt_t bsize, Float_t **data);
        ~TKDInterpolator();
 
-               TKDTreeIF* GetHelper() {return fKDhelper;}
+               Double_t   Eval(const Double_t *point, Double_t &result, Double_t &error);
        inline Bool_t     GetCOGPoint(Int_t node, Float_t *coord, Float_t &val, Float_t &error) const ;
        inline Bool_t     GetDataPoint(Int_t n, Float_t *p) const;
-               Double_t   Eval(const Double_t *point, Int_t npoints = 12);
-               void       DrawNodes(UInt_t ax1 = 0, UInt_t ax2 = 1, Int_t depth = -1);
+               TKDTreeIF* GetHelper() {return fKDhelper;}
+       inline Bool_t     GetIntInterpolation() const {return fStatus&1;}
+       inline Bool_t     GetSetStore() const {return fStatus&2;}
+       inline Bool_t     GetUseWeights() const {return fStatus&4;}
+                                       
+                                       void       DrawNodes(UInt_t ax1 = 0, UInt_t ax2 = 1, Int_t depth = -1);
                void       DrawNode(Int_t tnode, UInt_t ax1 = 0, UInt_t ax2 = 1);
-
+               void       GetStatus();
+               void       SetIntInterpolation(const Bool_t on = kTRUE);
+               void       SetSetStore(const Bool_t on = kTRUE);
+               void       SetUseWeights(const Bool_t on = kTRUE);
+       
 private:
        TKDInterpolator(const TKDInterpolator &);
        TKDInterpolator& operator=(const TKDInterpolator &);    
                void       Build();
-       
+               Double_t   CookPDF(const Double_t *point, const Int_t node, Double_t &result, Double_t &error);
+                                       
 protected:
        Int_t     fNTNodes;        //Number of evaluation data points
-       Float_t   **fRefPoints;    //[fNDim][fNTNodes]
+       Float_t   **fRefPoints;    //[fNDim]
        Float_t   *fRefValues;     //[fNTNodes]
+       TMatrixD  *fCov;           //[fNTNodes] cov matrix array for nodes
+       TVectorD  *fPar;           //[fNTNodes] parameters array for nodes
+       Bool_t    *fPDFstatus;     //[fNTNodes] status bit for node's PDF
 
 private:
-       Int_t                   fDepth;        //! depth of the KD Tree structure used
-       Double_t        *fTmpPoint;    //! temporary storage for one data point
-       TKDTreeIF *fKDhelper;    //! kNN finder
-       TLinearFitter *fFitter;  //! linear fitter      
+       UChar_t   fStatus;         // status of the interpolator
+       Int_t     fLambda;         // number of parameters in polynom
+       Int_t                   fDepth;          //! depth of the KD Tree structure used
+       Double_t        *fBuffer;        //! working space [2*fLambda]
+       TKDTreeIF *fKDhelper;      //! kNN finder
+       TLinearFitter *fFitter;    //! linear fitter    
 
        ClassDef(TKDInterpolator, 1)   // data interpolator based on KD tree
 };
@@ -48,7 +75,7 @@ Bool_t        TKDInterpolator::GetCOGPoint(Int_t node, Float_t *coord, Float_t &val, Fl
 
        for(int idim=0; idim<fNDim; idim++) coord[idim] = fRefPoints[idim][node];
        val   = fRefValues[node];
-       error = fRefValues[node]; // to be implemented
+       error = fRefValues[node]/TMath::Sqrt(fBucketSize);
        return kTRUE;
 }
 
index e679ed3..84ebccc 100644 (file)
@@ -5,6 +5,7 @@
 
 #ifndef R__ALPHA
 templateClassImp(TKDTree)
+templateClassImp(TKDNode)
 #endif
 //////////////////////////////////////////////////////////////////////
 // Memory setup of protected data members:
@@ -36,13 +37,18 @@ TKDTree<Index, Value>::TKDTree() :
        ,kDataOwner(kFALSE)
        ,fNnodes(0)
        ,fNDim(0)
+       ,fNDimm(0)
        ,fNpoints(0)
        ,fBucketSize(0)
-       ,fData(0x0)
+       ,fNodes(0x0)
        ,fRange(0x0)
+       ,fData(0x0)
        ,fBoundaries(0x0)
-       ,fNodes(0x0)
+       ,fkNNdim(0)
        ,fkNN(0x0)
+       ,fkNNdist(0x0)
+       ,fDistBuffer(0x0)
+       ,fIndBuffer(0x0)
        ,fIndPoints(0x0)
        ,fRowT0(0)
        ,fCrossNode(0)
@@ -59,20 +65,25 @@ TKDTree<Index, Value>::TKDTree(Index npoints, Index ndim, UInt_t bsize, Value **
        ,kDataOwner(kFALSE)
        ,fNnodes(0)
        ,fNDim(ndim)
+       ,fNDimm(2*ndim)
        ,fNpoints(npoints)
        ,fBucketSize(bsize)
-       ,fData(data) //Columnwise!!!!!
+       ,fNodes(0x0)
        ,fRange(0x0)
+       ,fData(data) //Columnwise!!!!!
        ,fBoundaries(0x0)
-       ,fNodes(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.
-       
+// Allocate data by hand. See TKDTree(TTree*, const Char_t*) for an automatic method.
+
        Build();
 }
 
@@ -80,6 +91,9 @@ TKDTree<Index, Value>::TKDTree(Index npoints, Index ndim, UInt_t bsize, Value **
 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 (fNodes) delete [] fNodes;
        if (fIndPoints) delete [] fIndPoints;
@@ -124,7 +138,7 @@ void TKDTree<Index, Value>::Build(){
        fRange = new Value[2*fNDim];
        fIndPoints= new Index[fNpoints];
        for (Index i=0; i<fNpoints; i++) fIndPoints[i] = i;
-       fNodes  = new TKDNode[fNnodes];
+       fNodes  = new TKDNode<Index, Value>[fNnodes];
        //
        fCrossNode = (1<<(fRowT0+1))-1;
        if (fCrossNode<fNnodes) fCrossNode = 2*fCrossNode+1;
@@ -166,7 +180,7 @@ void TKDTree<Index, Value>::Build(){
                Int_t crow     = rowStack[currentIndex];
                Int_t cpos     = posStack[currentIndex];
                Int_t cnode    = nodeStack[currentIndex];               
-               TKDNode * node = &(fNodes[cnode]);
+               TKDNode<Index, Value> * node = &(fNodes[cnode]);
                //printf("currentIndex %d npoints %d node %d\n", currentIndex, npoints, cnode);
                //
                // divide points
@@ -235,7 +249,7 @@ void TKDTree<Index, Value>::Build(){
 
 //_________________________________________________________________
 template <typename  Index, typename Value>
-Bool_t TKDTree<Index, Value>::FindNearestNeighbors(const Value *point, const Int_t k, Index *&in, Value &d)
+Int_t  TKDTree<Index, Value>::FindNearestNeighbors(const Value *point, const Int_t k, Index *&in, Value *&d)
 {
 // Find "k" nearest neighbors to "point".
 //
@@ -252,17 +266,31 @@ Bool_t    TKDTree<Index, Value>::FindNearestNeighbors(const Value *point, const Int
                return kFALSE;
        }
 
+       UInt_t debug = 0;
+       Int_t nCalculations = 0;
+       
        // 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. 
        Index nBounds = 0;
        if(!fBoundaries) MakeBoundaries();
@@ -271,9 +299,10 @@ Bool_t     TKDTree<Index, Value>::FindNearestNeighbors(const Value *point, const Int
                if(bounds[2*idim] == fRange[2*idim]) nBounds++;
                if(bounds[2*idim+1] == fRange[2*idim+1]) nBounds++;
        }
+       if(debug>=1) printf("Calculate boundaries [nBounds = %d]\n", nBounds);
        
        // traverse tree
-       TKDNode *node;
+       TKDNode<Index, Value> *node;
        Int_t nodeStack[128], nodeIn[128];
        Index currentIndex = 0;
        nodeStack[0]   = inode;
@@ -282,14 +311,15 @@ Bool_t    TKDTree<Index, Value>::FindNearestNeighbors(const Value *point, const Int
                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
                node = &fNodes[tnode/2 + (tnode%2) - 1];
-               if((TMath::Abs(point[node->fAxis] - node->fValue) > fkNN_dist[k-1])
+               if((TMath::Abs(point[node->fAxis] - node->fValue) > fkNNdist[k-1])
                                                                                        &&
                        ((point[node->fAxis] > node->fValue && tnode%2) ||
                         (point[node->fAxis] < node->fValue && !tnode%2))) {
-                       //printf("\tREMOVE NODE\n");
+                       //printf("\tREMOVE NODE %d\n", tnode/2 + (tnode%2) - 1);
 
                        // mark bound
                        nBounds++;
@@ -299,36 +329,47 @@ Bool_t    TKDTree<Index, Value>::FindNearestNeighbors(const Value *point, const Int
                }
                
                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 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
-                               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]);
+                                       }
+                               }                               
                        }
                }
                
@@ -337,7 +378,7 @@ Bool_t      TKDTree<Index, Value>::FindNearestNeighbors(const Value *point, const Int
                if(parent_node >= 0 && entry != parent_node){
                        // 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[node->fAxis] - node->fValue) > fkNNdist[k-1]){
                                // mark bound
                                nBounds++;
                                // end all recursions
@@ -346,14 +387,14 @@ Bool_t    TKDTree<Index, Value>::FindNearestNeighbors(const Value *point, const Int
                        currentIndex++;
                        nodeStack[currentIndex]=tnode/2 + (tnode%2) - 1;
                        nodeIn[currentIndex]=tnode;
-                       //printf("\tregister %d\n", nodeStack[currentIndex]);
+                       if(debug>=2) printf("\tregister %d\n", nodeStack[currentIndex]);
                }
-
+               
                // register children nodes
                Int_t child_node;
                Bool_t kAllow[] = {kTRUE, kTRUE};
                node = &fNodes[tnode];
-               if(TMath::Abs(point[node->fAxis] - node->fValue) > fkNN_dist[k-1]){
+               if(TMath::Abs(point[node->fAxis] - node->fValue) > fkNNdist[k-1]){
                        if(point[node->fAxis] > node->fValue) kAllow[0] = kFALSE;
                        else kAllow[1] = kFALSE;
                }
@@ -364,21 +405,15 @@ Bool_t    TKDTree<Index, Value>::FindNearestNeighbors(const Value *point, const Int
                                currentIndex++;
                                nodeStack[currentIndex] = child_node;
                                nodeIn[currentIndex]=tnode;
-                               //printf("\tregister %d\n", nodeStack[currentIndex]);
+                               if(debug>=2) printf("\tregister %d\n", nodeStack[currentIndex]);
                        }
-               }
+               }               
        }
        // 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 nCalculations; //kTRUE;
 }
 
 
@@ -398,7 +433,7 @@ Index TKDTree<Index, Value>::FindNode(const Value * point){
                currentIndex--;
                if (IsTerminal(inode)) return inode;
                
-               TKDNode & node = fNodes[inode];
+               TKDNode<Index, Value> & node = fNodes[inode];
                if (point[node.fAxis]<=node.fValue){
                        currentIndex++;
                        stackNode[currentIndex]=(inode*2)+1;
@@ -446,7 +481,7 @@ void TKDTree<Index, Value>::FindPoint(Value * point, Index &index, Int_t &iter){
                        continue;
                }
                
-               TKDNode & node = fNodes[inode];
+               TKDNode<Index, Value> & node = fNodes[inode];
                if (point[node.fAxis]<=node.fValue){
                        currentIndex++;
                        stackNode[currentIndex]=(inode*2)+1;
@@ -485,7 +520,7 @@ void TKDTree<Index, Value>::FindInRangeA(Value * point, Value * delta, Index *re
     currentIndex--;
     if (!IsTerminal(inode)){
       // not terminal
-      TKDNode * node = &(fNodes[inode]);
+      TKDNode<Index, Value> * node = &(fNodes[inode]);
       if (point[node->fAxis] - delta[node->fAxis] <  node->fValue) {
        currentIndex++; 
        stackNode[currentIndex]= (inode*2)+1;
@@ -573,7 +608,7 @@ void TKDTree<Index, Value>::FindInRangeB(Value * point, Value * delta, Index *re
       }
     }
     //
-    TKDNode * node = &(fNodes[inode]);
+    TKDNode<Index, Value> * node = &(fNodes[inode]);
     if (point[node->fAxis] - delta[node->fAxis] <  node->fValue) {
       currentIndex++; 
       stackNode[currentIndex]= (inode<<1)+1;
@@ -704,28 +739,28 @@ void TKDTree<Index, Value>::MakeBoundaries(Value *range)
        // total number of nodes including terminal nodes
        Int_t totNodes = fNnodes + GetNTerminalNodes();
        fBoundaries = new Value[totNodes*2*fNDim];
-       printf("Allocate %d nodes\n", totNodes);
+       Info("MakeBoundaries(Value*)", Form("Allocate boundaries for %d nodes", totNodes));
+       
        // loop
        Int_t child_index;
        for(int inode=fNnodes-1; inode>=0; inode--){
-               //printf("\tAllocate node %d\n", inode);
                memcpy(&fBoundaries[inode*2*fNDim], fRange, 2*fNDim*sizeof(Value));
                                
                // check left child node
                child_index = 2*inode+1;
-               if(IsTerminal(child_index)) CookBoundariesTerminal(inode, kTRUE);
+               if(IsTerminal(child_index)) CookBoundaries(inode, kTRUE);
                for(int idim=0; idim<fNDim; idim++) fBoundaries[2*fNDim*inode+2*idim] = fBoundaries[2*fNDim*child_index+2*idim];
                
                // check right child node
                child_index = 2*inode+2;
-               if(IsTerminal(child_index)) CookBoundariesTerminal(inode, kFALSE);
+               if(IsTerminal(child_index)) CookBoundaries(inode, kFALSE);
                for(int idim=0; idim<fNDim; idim++) fBoundaries[2*fNDim*inode+2*idim+1] = fBoundaries[2*fNDim*child_index+2*idim+1];
        }
 }
 
 //_________________________________________________________________
 template <typename Index, typename Value>
-void TKDTree<Index, Value>::CookBoundariesTerminal(Int_t parent_node, Bool_t LEFT)
+void TKDTree<Index, Value>::CookBoundaries(Int_t parent_node, Bool_t LEFT)
 {
        // define index of this terminal node
        Int_t index = 2 * parent_node + (LEFT ? 1 : 2);
@@ -738,7 +773,7 @@ void TKDTree<Index, Value>::CookBoundariesTerminal(Int_t parent_node, Bool_t LEF
        Int_t nvals = 0;
                
        // recurse parent nodes
-       TKDNode *node = 0x0;
+       TKDNode<Index, Value> *node = 0x0;
        while(parent_node >= 0 && nvals < 2 * fNDim){
                node = &(fNodes[parent_node]);
                if(LEFT){
index 3b48910..fb9866f 100644 (file)
@@ -6,6 +6,15 @@
 #endif
 
 #include "TMath.h"
+
+template <typename Index, typename Value>
+struct TKDNode{
+       Char_t fAxis;
+       Value  fValue;
+       
+       ClassDef(TKDNode, 1) // TKDTree node structure
+};
+
 template <typename Index, typename Value> class TKDTree : public TObject
 {
 public:
@@ -18,66 +27,65 @@ public:
        ~TKDTree();
        
        // getters
-                                       inline  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];
-                                       }
-                                       inline Char_t           GetNodeAxis(Int_t id) const {return (id < 0 || id >= fNnodes) ? 0 : fNodes[id].fAxis;}
-                                       inline Value            GetNodeValue(Int_t id) const {return (id < 0 || id >= fNnodes) ? 0 : fNodes[id].fValue;}
-                                       inline Int_t            GetNNodes() const {return fNnodes;}
-                                       inline Int_t            GetNTerminalNodes() const {return fNpoints/fBucketSize + ((fNpoints%fBucketSize)?1:0);}
-                                       inline Value*           GetBoundaries() const {return fBoundaries;}
-                                       inline Value*           GetBoundary(const Int_t node) const {return &fBoundaries[node*2*fNDim];}
+       inline  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];
+       }
+       inline Char_t           GetNodeAxis(Int_t id) const {return (id < 0 || id >= fNnodes) ? 0 : fNodes[id].fAxis;}
+       inline Value            GetNodeValue(Int_t id) const {return (id < 0 || id >= fNnodes) ? 0 : fNodes[id].fValue;}
+       inline Int_t            GetNNodes() const {return fNnodes;}
+       inline Int_t            GetNTerminalNodes() const {return fNpoints/fBucketSize + ((fNpoints%fBucketSize)?1:0);}
+       inline Value*           GetBoundaries();
+       inline Value*           GetBoundary(const Int_t node);
        static                                  Int_t           GetIndex(Int_t row, Int_t collumn){return collumn+(1<<row);}
        static  inline  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);
-                                       inline  void            FindBNodeA(Value * point, Value * delta, Int_t &inode);
-       //
-                                       inline  Bool_t  IsTerminal(Index inode){return (inode>=fNnodes);}
-       //
-                                                                       Value           KOrdStat(Index ntotal, Value *a, Index k, Index *index);
-                                                                       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);
+       Int_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);
+       inline  void    FindBNodeA(Value * point, Value * delta, Int_t &inode);
+       inline  Bool_t  IsTerminal(Index inode){return (inode>=fNnodes);}
+       Value           KOrdStat(Index ntotal, Value *a, Index k, Index *index);
+       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);
 
 protected:
-                                                                       void            Build();  // build tree
+       void            Build();  // build tree
                                                                        
 private:
-                                                                       TKDTree(const TKDTree &); // not implemented
-                                                                       TKDTree<Index, Value>& operator=(const TKDTree<Index, Value>&); // not implemented
-                                                                       void            CookBoundariesTerminal(Int_t parent_node, Bool_t left);
+       TKDTree(const TKDTree &); // not implemented
+       TKDTree<Index, Value>& operator=(const TKDTree<Index, Value>&); // not implemented
+       void            CookBoundaries(Int_t parent_node, Bool_t left);
 
-public:
-       struct TKDNode{
-               Char_t fAxis;
-               Value  fValue;
-       };
 
 protected:
-       Bool_t  kDataOwner;  // Toggle ownership of the data
+       Bool_t  kDataOwner;  //! Toggle ownership of the data
        Int_t   fNnodes;     // size of node array
-       Index   fNDim;       // number of variables
+       Index   fNDim;       // number of dimensions
+       Index   fNDimm;      // dummy 2*fNDim
        Index   fNpoints;    // number of multidimensional points
        Index   fBucketSize; // limit statistic for nodes 
-       Value   **fData;     //!
-       Value           *fRange;     //! range of data for each dimension     
+       TKDNode<Index, Value> *fNodes;     //[fNnodes] nodes descriptors array
+       Value           *fRange;     //[fNDimm] range of data for each dimension
+       Value   **fData;     //! data points
        Value           *fBoundaries;//! nodes boundaries - check class doc
-       TKDNode *fNodes;
-       Index   *fkNN;       //! k nearest neighbors
 
+// 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
+       Int_t   fRowT0;      //! smallest terminal row
+       Int_t   fCrossNode;  //! cross node
+       Int_t   fOffset;     //! offset in fIndPoints
 
        ClassDef(TKDTree, 1)  // KD tree
 };
@@ -87,17 +95,34 @@ typedef TKDTree<Int_t, Double_t> TKDTreeID;
 typedef TKDTree<Int_t, Float_t> TKDTreeIF;
 
 //_________________________________________________________________
-template <typename  Index, typename Value> void TKDTree<Index, Value>::FindBNodeA(Value *point, Value *delta, Int_t &inode){
+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;){
-    TKDNode & node = fNodes[inode];
+    TKDNode<Index, Value> & node = fNodes[inode];
     if (TMath::Abs(point[node.fAxis] - node.fValue)<delta[node.fAxis]) break;
     inode = (point[node.fAxis] < node.fValue) ? (inode*2)+1: (inode*2)+2;    
   }
 }
 
+//_________________________________________________________________
+template <typename  Index, typename Value>
+Value* TKDTree<Index, Value>::GetBoundaries()
+{
+       if(!fBoundaries) MakeBoundaries();
+       return fBoundaries;
+}
+
+//_________________________________________________________________
+template <typename  Index, typename Value>
+Value* TKDTree<Index, Value>::GetBoundary(const Int_t node)
+{
+       if(!fBoundaries) MakeBoundaries();
+       return &fBoundaries[node*2*fNDim];
+}
+
 #endif