#include "TRandom.h"
#include "TROOT.h"
-
-
ClassImp(TKDInterpolator)
+ClassImp(TKDInterpolator::TKDNodeInfo)
/////////////////////////////////////////////////////////////////////
// Memory setup of protected data memebers
//
// 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 )|
/////////////////////////////////////////////////////////////////////
//_________________________________________________________________
,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)
{
,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)
{
//_________________________________________________________________
-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)
{
//
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();
//_________________________________________________________________
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] ;
// - 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];
}
}
+//__________________________________________________________________
+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)
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.;
+}
+
#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
};
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;
}
#ifndef R__ALPHA
templateClassImp(TKDTree)
+templateClassImp(TKDNode)
#endif
//////////////////////////////////////////////////////////////////////
// Memory setup of protected data members:
,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)
,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();
}
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;
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;
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
//_________________________________________________________________
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".
//
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();
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;
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++;
}
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]);
+ }
+ }
}
}
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
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;
}
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;
}
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;
continue;
}
- TKDNode & node = fNodes[inode];
+ TKDNode<Index, Value> & node = fNodes[inode];
if (point[node.fAxis]<=node.fValue){
currentIndex++;
stackNode[currentIndex]=(inode*2)+1;
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;
}
}
//
- 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;
// 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);
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){
#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:
~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
};
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
git://git.uio.no / u / mrichter / AliRoot.git / commitdiff