rename Interpolator to PDF, add new class TKDInterpolator and
authorabercuci <abercuci@f7af4fe6-9843-0410-8265-dc069ae4e863>
Thu, 25 Oct 2007 14:04:29 +0000 (14:04 +0000)
committerabercuci <abercuci@f7af4fe6-9843-0410-8265-dc069ae4e863>
Thu, 25 Oct 2007 14:04:29 +0000 (14:04 +0000)
modify inheritance tree

14 files changed:
STAT/StatLinkDef.h
STAT/TKDInterpolator.cxx
STAT/TKDInterpolator.h
STAT/TKDInterpolatorBase.cxx [new file with mode: 0644]
STAT/TKDInterpolatorBase.h [new file with mode: 0644]
STAT/TKDNodeInfo.cxx [new file with mode: 0644]
STAT/TKDNodeInfo.h [new file with mode: 0644]
STAT/TKDPDF.cxx [new file with mode: 0644]
STAT/TKDPDF.h [new file with mode: 0644]
STAT/TKDSpline.cxx
STAT/TKDSpline.h
STAT/TKDTree.cxx
STAT/TKDTree.h
STAT/libSTAT.pkg

index dcbb258..3fa5253 100644 (file)
 #pragma link C++ typedef TKDTreeIF;
 
 
+#pragma link C++ class TKDInterpolatorBase+;
+#pragma link C++ class TKDNodeInfo+;
+#pragma link C++ class TKDPDF+;
 #pragma link C++ class TKDInterpolator+;
-#pragma link C++ class TKDInterpolator::TKDNodeInfo+;
 #pragma link C++ class TKDSpline+;
 
 #pragma link C++ class AliTMinuitToolkit+;
index 1d64336..b9d0c1f 100644 (file)
 #include "TKDInterpolator.h"
+#include "TKDNodeInfo.h"
 
-#include "TLinearFitter.h"
-#include "TTree.h"
-#include "TH2.h"
-#include "TObjArray.h"
-#include "TObjString.h"
-#include "TBox.h"
-#include "TGraph.h"
-#include "TMarker.h"
-#include "TVectorD.h"
-#include "TMatrixD.h"
+#include "TClonesArray.h"
 
 ClassImp(TKDInterpolator)
-ClassImp(TKDInterpolator::TKDNodeInfo)
 
-/////////////////////////////////////////////////////////////////////
-// Memory setup of protected data members
-// fRefPoints : evaluation point of PDF for each terminal node of underlying KD Tree.
-// | 1st terminal node (fNDim point coordinates) | 2nd terminal node (fNDim point coordinates) | ...
-//
-// 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 )|
-/////////////////////////////////////////////////////////////////////
 
-//_________________________________________________________________
-TKDInterpolator::TKDNodeInfo::TKDNodeInfo(const Int_t dim): 
-       fNDim(dim)
-       ,fRefPoint(0x0)
-       ,fRefValue(0.)
-       ,fCov(0x0)
-       ,fPar(0x0)
-{
-       if(fNDim) Build(dim);
-}
-
-//_________________________________________________________________
-TKDInterpolator::TKDNodeInfo::~TKDNodeInfo()
-{
-       if(fRefPoint) delete [] fRefPoint;
-       if(fCov){
-               delete fPar;
-               delete fCov;
-       }
-}
-
-//_________________________________________________________________
-void TKDInterpolator::TKDNodeInfo::Build(const Int_t dim)
-{
-// Allocate/Reallocate space for this node.
-
-       if(!dim) return;
-
-       fNDim = dim;
-       Int_t lambda = Int_t(1 + fNDim + .5*fNDim*(fNDim+1));
-       if(fRefPoint) delete [] fRefPoint;
-       fRefPoint = new Float_t[fNDim];
-       if(fCov){
-               fCov->ResizeTo(lambda, lambda);
-               fPar->ResizeTo(lambda);
-       }
-       return;
-}
 
 //_________________________________________________________________
-void TKDInterpolator::TKDNodeInfo::Store(const TVectorD &par, const TMatrixD &cov)
-{
-       if(!fCov){
-               fCov = new TMatrixD(cov.GetNrows(), cov.GetNrows());
-               fPar = new TVectorD(par.GetNrows());
-       }
-       (*fPar) = par;
-       (*fCov) = cov;
-}
-
-//_________________________________________________________________
-Double_t TKDInterpolator::TKDNodeInfo::CookPDF(const Double_t *point, Double_t &result, Double_t &error)
-{
-// Recalculate the PDF for one node from the results of interpolation (parameters and covariance matrix)
-
-       if(fNDim>10) return 0.; // support only up to 10 dimensions
-       
-       Int_t lambda = 1 + fNDim + (fNDim*(fNDim+1)>>1);
-       Double_t fdfdp[66];
-       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<lambda; i++){
-               result += fdfdp[i]*(*fPar)(i);
-               for(int j=0; j<lambda; j++) error += fdfdp[i]*fdfdp[j]*(*fCov)(i,j);
-       }       
-       error = TMath::Sqrt(error);
-       
-       //printf("TKDNodeInfo::CookPDF() : %6.3f +- %6.3f\n", result, error);
-
-       return 0.;
-}
-
-
-//_________________________________________________________________
-TKDInterpolator::TKDInterpolator() : TKDTreeIF()
-       ,fNTNodes(0)
-       ,fTNodes(0x0)
-       ,fStatus(4)
-       ,fLambda(0)
-       ,fDepth(-1)
-       ,fAlpha(.5)
-       ,fRefPoints(0x0)
-       ,fBuffer(0x0)
-       ,fKDhelper(0x0)
-       ,fFitter(0x0)
+TKDInterpolator::TKDInterpolator() :
+       TKDInterpolatorBase()
 {
 // Default constructor. To be used with care since in this case building
 // of data structure is completly left to the user responsability.
 }
 
 //_________________________________________________________________
-TKDInterpolator::TKDInterpolator(Int_t npoints, Int_t ndim, UInt_t bsize, Float_t **data) : TKDTreeIF(npoints, ndim, bsize, data)
-       ,fNTNodes(GetNTNodes())
-       ,fTNodes(0x0)
-       ,fStatus(4)
-       ,fLambda(0)
-       ,fDepth(-1)
-       ,fAlpha(.5)
-       ,fRefPoints(0x0)
-       ,fBuffer(0x0)
-       ,fKDhelper(0x0)
-       ,fFitter(0x0)
+TKDInterpolator::TKDInterpolator(Int_t ndim, Int_t npoints) :
+       TKDInterpolatorBase(ndim)
 {
 // Wrapper constructor for the TKDTree.
 
-       Build();
+       if(npoints) TKDInterpolatorBase::Build(npoints);
 }
 
 
 //_________________________________________________________________
-TKDInterpolator::TKDInterpolator(TTree *t, const Char_t *var, const Char_t *cut, UInt_t bsize, Long64_t nentries, Long64_t firstentry) : TKDTreeIF()
-       ,fNTNodes(0)
-       ,fTNodes(0x0)
-       ,fStatus(4)
-       ,fLambda(0)
-       ,fDepth(-1)
-       ,fAlpha(.5)
-       ,fRefPoints(0x0)
-       ,fBuffer(0x0)
-       ,fKDhelper(0x0)
-       ,fFitter(0x0)
-{
-// Alocate data from a tree. The variables which have to be analysed are
-// defined in the "var" parameter as a colon separated list. The format should
-// be identical to that used by TTree::Draw().
-//
-// 
-
-       TObjArray *vars = TString(var).Tokenize(":");
-       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", 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));
-                       fData = new Float_t*[fNDim];
-                       for(int idim=0; idim<fNDim; idim++) fData[idim] = new Float_t[fNpoints];
-                       kDataOwner = kTRUE;
-               }
-               v = t->GetV1();
-               for(int ip=0; ip<fNpoints; ip++) fData[idim][ip] = (Float_t)v[ip];
-       }
-       TKDTreeIF::Build();
-       Build();
-}
-
-//_________________________________________________________________
 TKDInterpolator::~TKDInterpolator()
 {
-       if(fFitter) delete fFitter;
-       if(fKDhelper) delete fKDhelper;
-       if(fBuffer) delete [] fBuffer;
-       
-       if(fRefPoints){
-               for(int idim=0; idim<fNDim; idim++) delete [] fRefPoints[idim] ;
-               delete [] fRefPoints;
-       }
-       if(fTNodes) delete [] fTNodes;
 }
 
 //_________________________________________________________________
-void TKDInterpolator::Build()
+void TKDInterpolator::AddNode(const TKDNodeInfo &node)
 {
-// Fill interpolator's data array i.e.
-//  - estimation points 
-//  - corresponding PDF values
-
-       fNTNodes = TKDTreeIF::GetNTNodes();
-       if(!fBoundaries) MakeBoundaries();
-       fLambda = 1 + fNDim + (fNDim*(fNDim+1)>>1);
-       //printf("after MakeBoundaries() %d\n", memory());
-       
-       // allocate memory for data
-       fTNodes = new TKDNodeInfo[fNTNodes];
-       for(int in=0; in<fNTNodes; in++) fTNodes[in].Build(fNDim);
-       //printf("after BuildNodes() %d\n", memory());
-
-       Float_t *bounds = 0x0;
-       Int_t *indexPoints;
-       for(int inode=0, tnode = fNnodes; inode<fNTNodes-1; inode++, tnode++){
-               fTNodes[inode].fRefValue =  Float_t(fBucketSize)/fNpoints;
-               bounds = GetBoundary(tnode);
-               for(int idim=0; idim<fNDim; idim++) fTNodes[inode].fRefValue /= (bounds[2*idim+1] - bounds[2*idim]);
-               
-               indexPoints = GetPointsIndexes(tnode);
-               // loop points in this terminal node
-               for(int idim=0; idim<fNDim; idim++){
-                       fTNodes[inode].fRefPoint[idim] = 0.;
-                       for(int ip = 0; ip<fBucketSize; ip++){
-/*                             printf("\t\tindex[%d] = %d %f\n", ip, indexPoints[ip], fData[idim][indexPoints[ip]]);*/
-                               fTNodes[inode].fRefPoint[idim] += fData[idim][indexPoints[ip]];
-                       }
-                       fTNodes[inode].fRefPoint[idim] /= fBucketSize;
-               }
-       }
-
-       // analyze last (incomplete) terminal node
-       Int_t counts = fNpoints%fBucketSize;
-       counts = counts ? counts : fBucketSize;
-       Int_t inode = fNTNodes - 1, tnode = inode + fNnodes;
-       fTNodes[inode].fRefValue =  Float_t(counts)/fNpoints;
-       bounds = GetBoundary(tnode);
-       for(int idim=0; idim<fNDim; idim++) fTNodes[inode].fRefValue /= (bounds[2*idim+1] - bounds[2*idim]);
-
-       // loop points in this terminal node
-       indexPoints = GetPointsIndexes(tnode);
-       for(int idim=0; idim<fNDim; idim++){
-               fTNodes[inode].fRefPoint[idim] = 0.;
-               for(int ip = 0; ip<counts; ip++) fTNodes[inode].fRefPoint[idim] += fData[idim][indexPoints[ip]];
-               fTNodes[inode].fRefPoint[idim] /= counts;
-       }
-}
-
-//__________________________________________________________________
-void TKDInterpolator::GetStatus()
-{
-// Prints the status of the interpolator
-
-       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");
-       return;
-       
-       printf("fNTNodes %d\n", fNTNodes);        //Number of evaluation data points
-       for(int i=0; i<fNTNodes; i++){
-               printf("%d ", i);
-               for(int idim=0; idim<fNDim; idim++) printf("%f ", fTNodes[i].fRefPoint[idim]);
-               printf("[%f] %s\n", fTNodes[i].fRefValue, fTNodes[i].fCov ? "true" : "false");
-               printf("Fit parameters : ");
-               if(!fTNodes[i].fPar){
-                       printf("Not defined.\n");
-                       continue;
-               }
-               for(int ip=0; ip<3; ip++) printf("p%d[%f] ", ip, (*fTNodes[i].fPar)(ip));
-               printf("\n");
-       }
-}
-
-//_________________________________________________________________
-Double_t TKDInterpolator::Eval(const Double_t *point, Double_t &result, Double_t &error, Bool_t force)
-{
-// 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((fStatus&1) && fTNodes[node].fCov && !force) return fTNodes[node].CookPDF(point, result, error);
-
-       // Allocate memory
-       if(!fBuffer) fBuffer = new Double_t[2*fLambda];
-       if(!fKDhelper){ 
-               fRefPoints = new Float_t*[fNDim];
-               for(int id=0; id<fNDim; id++){ 
-                       fRefPoints[id] = new Float_t[fNTNodes];
-                       for(int in=0; in<fNTNodes; in++) fRefPoints[id][in] = fTNodes[in].fRefPoint[id];
-               }
-//             for(int in=0; in<fNTNodes; in++){
-//                     printf("%3d ", in);
-//                     for(int id=0; id<fNDim; id++) printf("%6.3f ", fTNodes[in].fRefPoint[id]/*fRefPoints[id][in]*/);
-//                     printf("\n");
-//             }
-               fKDhelper = new TKDTreeIF(fNTNodes, fNDim, 30, fRefPoints);
-               fKDhelper->MakeBoundaries();
-       }
-       if(!fFitter) fFitter = new TLinearFitter(fLambda, Form("hyp%d", fLambda-1));
-       
-       // 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
-
-       Int_t *index,  // indexes of NN 
-             ipar,    // local looping variable
-                               npoints = Int_t((1.+fAlpha)*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;
-               }
-               // add points to fitter
-               fFitter->ClearPoints();
-               TKDNodeInfo *node = 0x0;
-               for(int in=0; in<npoints; in++){
-                       node = &fTNodes[index[in]];
-                       if(fStatus&1){ // INT
-                               Float_t *bounds = GetBoundary(FindNode(node->fRefPoint));
-                               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
-                               Float_t *p = node->fRefPoint;
-                               ipar = 0;
-                               for(int idim=0; idim<fNDim; idim++){
-                                       fBuffer[ipar++] = p[idim];
-                                       for(int jdim=idim; jdim<fNDim; jdim++) fBuffer[ipar++] = p[idim]*p[jdim];
-                               }
-                       }
-
-                       // 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, node->fRefValue, node->fRefValue*sig/w);
-               }
-               npoints += 4;
-       } while(fFitter->Eval());
-
-       // 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) fTNodes[node].Store(par, cov);
-               
-       // Build df/dpi|x values
-       Double_t *fdfdp = &fBuffer[fLambda];
-       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];
+       if(!fTNodes){
+               printf("W - TKDInterpolator::SetNode() : Node array not defined.\n");
+               return;
        }
 
-       // 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;
+       new((*fTNodes)[fNTNodes++]) TKDNodeInfo(node);
 }
 
 //_________________________________________________________________
-void TKDInterpolator::DrawNodes(UInt_t ax1, UInt_t ax2, Int_t depth)
+void TKDInterpolator::Build(Int_t npoints, Int_t ndim)
 {
-// Draw nodes structure projected on plane "ax1:ax2". The parameter
-// "depth" specifies the bucket size per node. If depth == -1 draw only
-// terminal nodes and evaluation points (default -1 i.e. bucket size per node equal bucket size specified by the user)
-//
-// Observation:
-// This function creates the nodes (TBox) array for the specified depth
-// but don't delete it. Abusing this function may cause memory leaks !
-
-
-       if(!fBoundaries) MakeBoundaries();
-
-       // Count nodes in specific view
-       Int_t nnodes = 0;
-       for(int inode = 0; inode <= 2*fNnodes; inode++){
-               if(depth == -1){
-                       if(!IsTerminal(inode)) continue;
-               } else if((inode+1) >> depth != 1) continue;
-               nnodes++;
-       }
-
-       //printf("depth %d nodes %d\n", depth, nnodes);
-       
-       TH2 *h2 = new TH2S("hNodes", "", 100, fRange[2*ax1], fRange[2*ax1+1], 100, fRange[2*ax2], fRange[2*ax2+1]);
-       h2->GetXaxis()->SetTitle(Form("x_{%d}", ax1));
-       h2->GetYaxis()->SetTitle(Form("x_{%d}", ax2));
-       h2->Draw();
-       
-       const Float_t kBorder = 0.;//1.E-4;
-       TBox *nodeArray = new TBox[nnodes], *node;
-       Float_t *bounds = 0x0;
-       nnodes = 0;
-       for(int inode = 0; inode <= 2*fNnodes; inode++){
-               if(depth == -1){
-                       if(!IsTerminal(inode)) continue;
-               } else if((inode+1) >> depth != 1) continue;
-
-               node = &nodeArray[nnodes++];
-               //node = new TBox(bounds[2*ax1]+border, bounds[2*ax2]+border, bounds[2*ax1+1]-border, bounds[2*ax2+1]-border);
-               node->SetFillStyle(3002);       
-               node->SetFillColor(50+inode/*Int_t(gRandom->Uniform()*50.)*/);
-               bounds = GetBoundary(inode);
-               node->DrawBox(bounds[2*ax1]+kBorder, bounds[2*ax2]+kBorder, bounds[2*ax1+1]-kBorder, bounds[2*ax2+1]-kBorder);
-       }
-       if(depth != -1) return;
-
-       // Draw reference points
-       TGraph *ref = new TGraph(fNTNodes);
-       ref->SetMarkerStyle(3);
-       ref->SetMarkerSize(.7);
-       ref->SetMarkerColor(2);
-       for(int inode = 0; inode < fNTNodes; inode++) ref->SetPoint(inode, fTNodes[inode].fRefPoint[ax1], fTNodes[inode].fRefPoint[ax2]);
-       ref->Draw("p");
-       return;
+       fNSize = ndim;
+       TKDInterpolatorBase::Build(npoints);
 }
 
 //_________________________________________________________________
-void TKDInterpolator::DrawNode(Int_t tnode, UInt_t ax1, UInt_t ax2)
-{
-// Draw node "node" and the data points within.
-//
-// Observation:
-// This function creates some graphical objects
-// but don't delete it. Abusing this function may cause memory leaks !
-
-       if(tnode < 0 || tnode >= fNTNodes){
-               Warning("DrawNode()", Form("Terminal node %d outside defined range.", tnode));
-               return;
+Int_t TKDInterpolator::GetNodeIndex(const Float_t *p)
+{
+/*     printf("TKDInterpolator::GetNodeIndex() ...\n");
+       printf("Looking for p[");
+       for(int i=0; i<fNSize; i++) printf("%f ", p[i]);
+       printf("] ...\n");*/
+       TKDNodeInfo *node;
+       for(int inode=0; inode<fNTNodes; inode++){
+               node = (TKDNodeInfo*)(*fTNodes)[inode];
+               //node->Print();
+               if(node->Has(p)) return inode;
        }
-
-       Int_t inode = tnode;
-       tnode += fNnodes;
-       // select zone of interest in the indexes array
-       Int_t *index = GetPointsIndexes(tnode);
-       Int_t nPoints = (tnode == 2*fNnodes) ? fNpoints%fBucketSize : fBucketSize;
-
-       // draw data points
-       TGraph *g = new TGraph(nPoints);
-       g->SetMarkerStyle(7);
-       for(int ip = 0; ip<nPoints; ip++) g->SetPoint(ip, fData[ax1][index[ip]], fData[ax2][index[ip]]);
-
-       // draw estimation point
-       TMarker *m=new TMarker(fTNodes[inode].fRefPoint[ax1], fTNodes[inode].fRefPoint[ax2], 20);
-       m->SetMarkerColor(2);
-       m->SetMarkerSize(1.7);
-       
-       // draw node contour
-       Float_t *bounds = GetBoundary(tnode);
-       TBox *n = new TBox(bounds[2*ax1], bounds[2*ax2], bounds[2*ax1+1], bounds[2*ax2+1]);
-       n->SetFillStyle(0);
-
-       g->Draw("ap");
-       m->Draw();
-       n->Draw();
-       
-       return;
-}
-
-
-//__________________________________________________________________
-void TKDInterpolator::SetInterpolationMethod(const Bool_t on)
-{
-// Set interpolation bit to "on".
-       
-       if(on) fStatus += fStatus&1 ? 0 : 1;
-       else fStatus += fStatus&1 ? -1 : 0;
+       return -1;
 }
 
 
 //_________________________________________________________________
-void TKDInterpolator::SetStore(const Bool_t on)
+Bool_t TKDInterpolator::SetNode(const Int_t inode, const TKDNodeInfo &ref)
 {
-// Set store bit to "on"
-       
-       if(on) fStatus += fStatus&2 ? 0 : 2;
-       else fStatus += fStatus&2 ? -2 : 0;
+       if(!fTNodes){
+               printf("W - TKDInterpolator::SetNode() : Node array not defined.\n");
+               return kFALSE;
+       }
+       if(inode >= fNTNodes){
+               printf("W - TKDInterpolator::SetNode() : Node array defined up to %d.\n", fNTNodes);
+               return kFALSE;
+       }
+       TKDNodeInfo *node = (TKDNodeInfo*)(*fTNodes)[inode];
+       (*node) = ref;
 }
 
-//_________________________________________________________________
-void TKDInterpolator::SetWeights(const Bool_t on)
-{
-// Set weights bit to "on"
-       
-       if(on) fStatus += fStatus&4 ? 0 : 4;
-       else fStatus += fStatus&4 ? -4 : 0;
-}
index b8913bf..94b4783 100644 (file)
 #ifndef ROOT_TKDInterpolator
 #define ROOT_TKDInterpolator
 
-#ifndef ROOT_TKDTree
-#include "TKDTree.h"
+#ifndef ROOT_TKDInterpolatorBase
+#include "TKDInterpolatorBase.h"
 #endif
 
-// Non parametric interpolation class based on local polinomial regression.
-// The class can be used to approximate PDF together with TKDTree or for
-// general regression when the data points are given.
-
-template <typename Value> class TVectorT;
-typedef struct TVectorT<Double_t> TVectorD;
-template <typename Value> class TMatrixT;
-typedef class TMatrixT<Double_t> TMatrixD;
-template <typename Index, typename Value> class TKDTree;
-typedef class TKDTree<Int_t, Float_t> TKDTreeIF;
-class TTree;
-class TLinearFitter;
-class TKDInterpolator : public TKDTreeIF
+class TKDInterpolator : public TKDInterpolatorBase
 {
 public:
        TKDInterpolator();
-       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(Int_t ndim, Int_t npoints=0);
        ~TKDInterpolator();
+       void       AddNode(const TKDNodeInfo &ref);
+       void       Build(Int_t npoints, Int_t ndim);
+       Int_t      GetNodeIndex(const Float_t *p);
+       Bool_t     SetNode(const Int_t i, const TKDNodeInfo &ref);
 
-               Double_t   Eval(const Double_t *point, Double_t &result, Double_t &error, Bool_t force = kFALSE);
-                                       Float_t    GetAlpha() const {return fAlpha;}
-       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;
-               TKDTreeIF* GetHelper() {return fKDhelper;}
-               Bool_t     GetInterpolationMethod() const {return fStatus&1;}
-               Int_t      GetNTNodes() const {return fNTNodes;}
-               Bool_t     GetStore() const {return fStatus&2;}
-               Bool_t     GetWeights() 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       SetAlpha(const Float_t a){if(a>0.) fAlpha = a;}
-               void       SetInterpolationMethod(const Bool_t on = kTRUE);
-               void       SetStore(const Bool_t on = kTRUE);
-               void       SetWeights(const Bool_t on = kTRUE);
-       
 private:
        TKDInterpolator(const TKDInterpolator &);
        TKDInterpolator& operator=(const TKDInterpolator &);    
-               void       Build();
-                                       
-public:
-       class TKDNodeInfo
-       {
-       public:
-               TKDNodeInfo(const Int_t ndim = 0);
-               virtual  ~TKDNodeInfo();
-               void      Build(const Int_t ndim);
-         Double_t  CookPDF(const Double_t *point, Double_t &result, Double_t &error);
-               void      Store(const TVectorD &par, const TMatrixD &cov);
-       
-               Int_t     fNDim;         // data dimension
-               Float_t   *fRefPoint;    //[fNDim] node's COG
-               Float_t   fRefValue;     // measured value for node 
-               TMatrixD  *fCov;         // interpolator covariance matrix
-               TVectorD  *fPar;         // interpolator parameters
-
-       private:
-               TKDNodeInfo(const TKDNodeInfo &);
-               TKDNodeInfo& operator=(const TKDNodeInfo &);
-               
-               ClassDef(TKDNodeInfo, 1) // node info for interpolator
-       };
-
-protected:
-       Int_t     fNTNodes;        //Number of evaluation data points
-       TKDNodeInfo *fTNodes;      //[fNTNodes] interpolation node
 
 private:
-       UChar_t   fStatus;         // status of the interpolator
-       UChar_t   fLambda;         // number of parameters in polynom
-       Short_t         fDepth;          //! depth of the KD Tree structure used
-       Float_t   fAlpha;          // alpha parameter
-       Float_t   **fRefPoints;    //! temporary storage of COG data
-       Double_t        *fBuffer;        //! working space [2*fLambda]
-       TKDTreeIF *fKDhelper;      //! kNN finder
-       TLinearFitter *fFitter;    //! linear fitter    
-
-       ClassDef(TKDInterpolator, 1)   // data interpolator based on KD tree
+       
+       ClassDef(TKDInterpolator, 1)   // LOWESS data interpolator
 };
 
-//__________________________________________________________________
-Bool_t TKDInterpolator::GetCOGPoint(Int_t node, Float_t *coord, Float_t &val, Float_t &err) const
-{
-       if(node < 0 || node > fNTNodes) return kFALSE;
-
-       for(int idim=0; idim<fNDim; idim++) coord[idim] = fTNodes[node].fRefPoint[idim];
-       val = fTNodes[node].fRefValue;
-       err = fTNodes[node].fRefValue/TMath::Sqrt(fBucketSize);
-       return kTRUE;
-}
-
-//__________________________________________________________________
-Bool_t TKDInterpolator::GetDataPoint(Int_t n, Float_t *p) const
-{
-       if(n < 0 || n >= fNpoints) return kFALSE;
-       if(!fData) return kFALSE;
-               
-       for(int i=0; i<fNDim; i++) p[i] = fData[i][n];
-       return kTRUE;
-}
-
 
 #endif
 
diff --git a/STAT/TKDInterpolatorBase.cxx b/STAT/TKDInterpolatorBase.cxx
new file mode 100644 (file)
index 0000000..7b1d61b
--- /dev/null
@@ -0,0 +1,307 @@
+#include "TKDInterpolatorBase.h"
+#include "TKDNodeInfo.h"
+#include "TKDTree.h"
+
+#include "TClonesArray.h"
+#include "TLinearFitter.h"
+#include "TTree.h"
+#include "TH2.h"
+#include "TObjArray.h"
+#include "TObjString.h"
+#include "TBox.h"
+#include "TGraph.h"
+#include "TMarker.h"
+#include "TMath.h"
+#include "TVectorD.h"
+#include "TMatrixD.h"
+
+ClassImp(TKDInterpolatorBase)
+
+/////////////////////////////////////////////////////////////////////
+// Memory setup of protected data members
+// fRefPoints : evaluation point of PDF for each terminal node of underlying KD Tree.
+// | 1st terminal node (fNDim point coordinates) | 2nd terminal node (fNDim point coordinates) | ...
+//
+// 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 )|
+/////////////////////////////////////////////////////////////////////
+
+
+//_________________________________________________________________
+TKDInterpolatorBase::TKDInterpolatorBase(const Int_t dim) :
+       fNSize(dim)
+       ,fNTNodes(0)
+       ,fTNodes(0x0)
+       ,fStatus(4)
+       ,fLambda(1 + dim + (dim*(dim+1)>>1))
+       ,fDepth(-1)
+       ,fAlpha(.5)
+       ,fRefPoints(0x0)
+       ,fBuffer(0x0)
+       ,fKDhelper(0x0)
+       ,fFitter(0x0)
+{
+// Default constructor. To be used with care since in this case building
+// of data structure is completly left to the user responsability.
+}
+
+//_________________________________________________________________
+void   TKDInterpolatorBase::Build(const Int_t n)
+{
+       // allocate memory for data
+
+       if(fTNodes) delete fTNodes;
+       fNTNodes = n;
+       fTNodes = new TClonesArray("TKDNodeInfo", fNTNodes);
+       for(int in=0; in<fNTNodes; in++) new ((*fTNodes)[in]) TKDNodeInfo(fNSize);
+}
+
+//_________________________________________________________________
+TKDInterpolatorBase::~TKDInterpolatorBase()
+{
+       if(fFitter) delete fFitter;
+       if(fKDhelper) delete fKDhelper;
+       if(fBuffer) delete [] fBuffer;
+       
+       if(fRefPoints){
+               for(int idim=0; idim<fNSize; idim++) delete [] fRefPoints[idim] ;
+               delete [] fRefPoints;
+       }
+       if(fTNodes) delete fTNodes;
+}
+
+
+//__________________________________________________________________
+Bool_t TKDInterpolatorBase::GetCOGPoint(Int_t inode, Float_t *&coord, Float_t &val, Float_t &err) const
+{
+       if(inode < 0 || inode > fNTNodes) return kFALSE;
+
+       TKDNodeInfo *node = (TKDNodeInfo*)(*fTNodes)[inode];
+       coord = &(node->fData[0]);
+       val = node->fVal[0];
+       err = node->fVal[1];
+       return kTRUE;
+}
+
+//_________________________________________________________________
+TKDNodeInfo* TKDInterpolatorBase::GetNodeInfo(const Int_t inode) const
+{
+       if(!fTNodes || inode >= fNTNodes) return 0x0;
+       return (TKDNodeInfo*)(*fTNodes)[inode];
+}
+
+
+//__________________________________________________________________
+void TKDInterpolatorBase::GetStatus()
+{
+// Prints the status of the interpolator
+
+       printf("Interpolator Status :\n");
+       printf("  Dim    : %d [%d]\n", fNSize, fLambda);
+       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("fNTNodes %d\n", fNTNodes);        //Number of evaluation data points
+       for(int i=0; i<fNTNodes; i++){
+               TKDNodeInfo *node = (TKDNodeInfo*)(*fTNodes)[i]; 
+               printf("%d ", i); node->Print();
+       }
+}
+
+//_________________________________________________________________
+Double_t TKDInterpolatorBase::Eval(const Double_t *point, Double_t &result, Double_t &error, Bool_t force)
+{
+// 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<fNSize; idim++) pointF[idim] = (Float_t)point[idim];
+       Int_t nodeIndex = GetNodeIndex(pointF);
+       if(nodeIndex<0){
+               result = 0.;
+               error = 1.E10;
+               return 0.;
+       }
+       TKDNodeInfo *node = (TKDNodeInfo*)(*fTNodes)[nodeIndex];
+       if((fStatus&1) && node->fCov && !force) return node->CookPDF(point, result, error);
+
+       // Allocate memory
+       if(!fBuffer) fBuffer = new Double_t[2*fLambda];
+       if(!fKDhelper){ 
+               fRefPoints = new Float_t*[fNSize];
+               for(int id=0; id<fNSize; id++){
+                       fRefPoints[id] = new Float_t[fNTNodes];
+                       for(int in=0; in<fNTNodes; in++) fRefPoints[id][in] = ((TKDNodeInfo*)(*fTNodes)[in])->fData[id];
+               }
+               fKDhelper = new TKDTreeIF(fNTNodes, fNSize, 30, fRefPoints);
+               fKDhelper->MakeBoundaries();
+       }
+       if(!fFitter) fFitter = new TLinearFitter(fLambda, Form("hyp%d", fLambda-1));
+       
+       // 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
+
+       Int_t *index,  // indexes of NN 
+             ipar,    // local looping variable
+                               npoints = Int_t((1.+fAlpha)*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
+
+       do{
+               // find nearest neighbors
+               for(int idim=0; idim<fNSize; 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<fNSize; idim++) printf("%f ", point[idim]);
+                       printf("\n");
+                       return -1;
+               }
+               // add points to fitter
+               fFitter->ClearPoints();
+               TKDNodeInfo *tnode = 0x0;
+               for(int in=0; in<npoints; in++){
+                       tnode = (TKDNodeInfo*)(*fTNodes)[index[in]];
+                       //tnode->Print();
+                       if(fStatus&1){ // INT
+                               Float_t *bounds = &(tnode->fData[fNSize]);
+                               ipar = 0;
+                               for(int idim=0; idim<fNSize; 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<fNSize; jdim++) fBuffer[ipar++] = (bounds[2*idim] + bounds[2*idim+1]) * (bounds[2*jdim] + bounds[2*jdim+1]) * .25;
+                               }
+                       } else { // COG
+                               Float_t *p = &(tnode->fData[0]);
+                               ipar = 0;
+                               for(int idim=0; idim<fNSize; idim++){
+                                       fBuffer[ipar++] = p[idim];
+                                       for(int jdim=idim; jdim<fNSize; jdim++) fBuffer[ipar++] = p[idim]*p[jdim];
+                               }
+                       }
+
+                       // 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.;
+                        
+//                     printf("x[");
+//                     for(int idim=0; idim<fLambda-1; idim++) printf("%f ", fBuffer[idim]);
+//                     printf("]  v[%f +- %f] (%f, %f)\n", tnode->fVal[0], tnode->fVal[1]/w, tnode->fVal[1], w);
+                       fFitter->AddPoint(fBuffer, tnode->fVal[0], tnode->fVal[1]/w);
+               }
+               npoints += 4;
+       } while(fFitter->Eval());
+
+       // 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) node->Store(par, cov);
+               
+       // Build df/dpi|x values
+       Double_t *fdfdp = &fBuffer[fLambda];
+       ipar = 0;
+       fdfdp[ipar++] = 1.;
+       for(int idim=0; idim<fNSize; idim++){
+               fdfdp[ipar++] = point[idim];
+               for(int jdim=idim; jdim<fNSize; jdim++) fdfdp[ipar++] = point[idim]*point[jdim];
+       }
+
+       // 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;
+}
+
+//_________________________________________________________________
+void TKDInterpolatorBase::DrawBins(UInt_t ax1, UInt_t ax2, Float_t ax1min, Float_t ax1max, Float_t ax2min, Float_t ax2max)
+{
+// Draw nodes structure projected on plane "ax1:ax2". The parameter
+// "depth" specifies the bucket size per node. If depth == -1 draw only
+// terminal nodes and evaluation points (default -1 i.e. bucket size per node equal bucket size specified by the user)
+//
+// Observation:
+// This function creates the nodes (TBox) array for the specified depth
+// but don't delete it. Abusing this function may cause memory leaks !
+
+
+       
+       TH2 *h2 = new TH2S("hNodes", "", 100, ax1min, ax1max, 100, ax2min, ax2max);
+       h2->GetXaxis()->SetTitle(Form("x_{%d}", ax1));
+       h2->GetYaxis()->SetTitle(Form("x_{%d}", ax2));
+       h2->Draw();
+       
+       const Float_t kBorder = 0.;//1.E-4;
+       TBox *boxArray = new TBox[fNTNodes], *box;
+       Float_t *bounds = 0x0;
+       for(int inode = 0; inode < fNTNodes; inode++){
+               box = &boxArray[inode];
+               box->SetFillStyle(3002);
+               box->SetFillColor(50+inode/*Int_t(gRandom->Uniform()*50.)*/);
+               
+               bounds = &(((TKDNodeInfo*)(*fTNodes)[inode])->fData[fNSize]);
+               box->DrawBox(bounds[2*ax1]+kBorder, bounds[2*ax2]+kBorder, bounds[2*ax1+1]-kBorder, bounds[2*ax2+1]-kBorder);
+       }
+
+       // Draw reference points
+       TGraph *ref = new TGraph(fNTNodes);
+       ref->SetMarkerStyle(3);
+       ref->SetMarkerSize(.7);
+       ref->SetMarkerColor(2);
+       for(int inode = 0; inode < fNTNodes; inode++){
+               TKDNodeInfo *node = (TKDNodeInfo*)(*fTNodes)[inode];
+               ref->SetPoint(inode, node->fData[ax1], node->fData[ax2]);
+       }
+       ref->Draw("p");
+       return;
+}
+
+//__________________________________________________________________
+void TKDInterpolatorBase::SetInterpolationMethod(const Bool_t on)
+{
+// Set interpolation bit to "on".
+       
+       if(on) fStatus += fStatus&1 ? 0 : 1;
+       else fStatus += fStatus&1 ? -1 : 0;
+}
+
+
+//_________________________________________________________________
+void TKDInterpolatorBase::SetStore(const Bool_t on)
+{
+// Set store bit to "on"
+       
+       if(on) fStatus += fStatus&2 ? 0 : 2;
+       else fStatus += fStatus&2 ? -2 : 0;
+}
+
+//_________________________________________________________________
+void TKDInterpolatorBase::SetWeights(const Bool_t on)
+{
+// Set weights bit to "on"
+       
+       if(on) fStatus += fStatus&4 ? 0 : 4;
+       else fStatus += fStatus&4 ? -4 : 0;
+}
diff --git a/STAT/TKDInterpolatorBase.h b/STAT/TKDInterpolatorBase.h
new file mode 100644 (file)
index 0000000..5b59170
--- /dev/null
@@ -0,0 +1,76 @@
+#ifndef ROOT_TKDInterpolatorBase
+#define ROOT_TKDInterpolatorBase
+
+#ifndef ROOT_Rtypes
+#include "Rtypes.h"
+#endif
+
+///////////////////////////////////////////////////////////////
+//
+// Base non parametric interpolation algorithm.
+// The class implements local polynomial regression (LOWESS).
+// The user will work with daughter classes which implements
+// particular data configurations.
+//
+///////////////////////////////////////////////////////////////
+
+template <typename Value> class TVectorT;
+typedef struct TVectorT<Double_t> TVectorD;
+template <typename Value> class TMatrixT;
+typedef class TMatrixT<Double_t> TMatrixD;
+template <typename Index, typename Value> class TKDTree;
+typedef class TKDTree<Int_t, Float_t> TKDTreeIF;
+class TLinearFitter;
+class TClonesArray;
+class TKDNodeInfo;
+class TKDInterpolatorBase
+{
+public:        
+       TKDInterpolatorBase(const Int_t size = 0);
+       virtual ~TKDInterpolatorBase();
+
+       Double_t   Eval(const Double_t *point, Double_t &result, Double_t &error, Bool_t force = kFALSE);
+       virtual Int_t GetNodeIndex(const Float_t *p) = 0;
+       Float_t    GetAlpha() const {return fAlpha;}
+       Bool_t     GetCOGPoint(Int_t node, Float_t *&coord, Float_t &val, Float_t &error) const;
+       Bool_t     GetInterpolationMethod() const {return fStatus&1;}
+       TKDNodeInfo* GetNodeInfo(const Int_t inode) const;
+       Int_t      GetNTNodes() const {return fNTNodes;}
+       void       GetStatus();
+       Bool_t     GetStore() const {return fStatus&2;}
+       Bool_t     GetWeights() const {return fStatus&4;}
+
+       void       DrawBins(UInt_t ax1 = 0, UInt_t ax2 = 1, Float_t ax1min=-1., Float_t ax1max=1., Float_t ax2min=-1., Float_t ax2max=1.);
+       void       SetAlpha(const Float_t a){if(a>0.) fAlpha = a;}
+       void       SetInterpolationMethod(const Bool_t on = kTRUE);
+       void       SetStore(const Bool_t on = kTRUE);
+       void       SetWeights(const Bool_t on = kTRUE);
+
+protected:
+       virtual void      Build(const Int_t nnodes);
+
+private:
+       TKDInterpolatorBase(const TKDInterpolatorBase &);
+       TKDInterpolatorBase& operator=(const TKDInterpolatorBase &);
+
+protected:
+       Int_t         fNSize;       // data dimension
+       Int_t         fNTNodes;     //Number of evaluation data points
+       TClonesArray  *fTNodes;     //interpolation nodes
+
+//private:
+       UChar_t       fStatus;      // status of the interpolator
+       UChar_t       fLambda;      // number of parameters in polynom
+       Short_t             fDepth;       //! depth of the KD Tree structure used
+       Float_t       fAlpha;       // alpha parameter
+       Float_t       **fRefPoints; //! temporary storage of COG data
+       Double_t            *fBuffer;     //! working space [2*fLambda]
+       TKDTree<Int_t, Float_t> *fKDhelper;      //! kNN finder
+       TLinearFitter *fFitter;     //! linear fitter   
+
+       ClassDef(TKDInterpolatorBase, 1)   // data interpolator based on KD tree
+};
+
+
+#endif
+
diff --git a/STAT/TKDNodeInfo.cxx b/STAT/TKDNodeInfo.cxx
new file mode 100644 (file)
index 0000000..a406952
--- /dev/null
@@ -0,0 +1,152 @@
+#include "TKDNodeInfo.h"
+
+#include "TVectorD.h"
+#include "TMatrixD.h"
+#include "TMath.h"
+
+ClassImp(TKDNodeInfo)
+
+
+//_________________________________________________________________
+TKDNodeInfo::TKDNodeInfo(const Int_t dim):
+       TObject()
+       ,fNDim(3*dim)
+       ,fData(0x0)
+       ,fCov(0x0)
+       ,fPar(0x0)
+{
+       fVal[0] = 0.; fVal[1] = 0.;
+       Build(dim);
+}
+
+//_________________________________________________________________
+TKDNodeInfo::TKDNodeInfo(const TKDNodeInfo &ref):
+       TObject((TObject&) ref)
+       ,fNDim(fNDim)
+       ,fData(0x0)
+       ,fCov(0x0)
+       ,fPar(0x0)
+{
+       Build(fNDim/3);
+
+       memcpy(fData, ref.fData, fNDim*sizeof(Float_t));
+       fVal[0] = ref.fVal[0];
+       fVal[1] = ref.fVal[1];
+       if(ref.fCov) (*fCov) = (*ref.fCov);
+       if(ref.fPar) (*fPar) = (*ref.fPar);
+}
+
+
+//_________________________________________________________________
+TKDNodeInfo::~TKDNodeInfo()
+{
+       if(fData) delete [] fData;
+       if(fCov){
+               delete fPar;
+               delete fCov;
+       }
+}
+
+//_________________________________________________________________
+TKDNodeInfo& TKDNodeInfo::operator=(const TKDNodeInfo & ref)
+{
+//     Info("operator==()", "...");
+       
+       Int_t ndim = fNDim/3;
+       if(fNDim != ref.fNDim){
+               fNDim = ref.fNDim;
+               Build(ndim);
+       }
+       memcpy(fData, ref.fData, fNDim*sizeof(Float_t));
+       fVal[0] = ref.fVal[0];
+       fVal[1] = ref.fVal[1];
+       if(ref.fCov) (*fCov) = (*ref.fCov);
+       if(ref.fPar) (*fPar) = (*ref.fPar);
+       
+       return *this;
+}
+
+//_________________________________________________________________
+void TKDNodeInfo::Build(const Int_t dim)
+{
+// Allocate/Reallocate space for this node.
+
+//     Info("Build()", "...");
+
+       if(!dim) return;
+       
+       Int_t lambda = Int_t(1 + dim + .5*dim*(dim+1));
+       if(fData) delete [] fData;
+       fData = new Float_t[fNDim];
+       if(fCov){
+               fCov->ResizeTo(lambda, lambda);
+               fPar->ResizeTo(lambda);
+       }
+       return;
+}
+
+//_________________________________________________________________
+void TKDNodeInfo::Print()
+{
+       Int_t dim = fNDim/3;
+       printf("x[");
+       for(int idim=0; idim<dim; idim++) printf("%f ", fData[idim]);
+       printf("] f = [%f +- %f]\n", fVal[0], fVal[1]);
+
+       Float_t *bounds = &fData[dim];
+       printf("range[");
+       for(int idim=0; idim<dim; idim++) printf("(%f %f) ", fData[2*idim], fData[2*idim+1]);
+       printf("]\n");
+       
+       printf("Fit parameters : ");
+       if(!fCov){
+               printf("Not defined.\n");
+               return;
+       }
+       
+       Int_t lambda = Int_t(1 + dim + .5*dim*(dim+1));
+       for(int ip=0; ip<3; ip++) printf("p%d[%f] ", ip, (*fPar)(ip));
+       printf("\n");
+}
+
+//_________________________________________________________________
+void TKDNodeInfo::Store(const TVectorD &par, const TMatrixD &cov)
+{
+       if(!fCov){
+               fCov = new TMatrixD(cov.GetNrows(), cov.GetNrows());
+               fPar = new TVectorD(par.GetNrows());
+       }
+       (*fPar) = par;
+       (*fCov) = cov;
+}
+
+//_________________________________________________________________
+Double_t TKDNodeInfo::CookPDF(const Double_t *point, Double_t &result, Double_t &error)
+{
+// Recalculate the PDF for one node from the results of interpolation (parameters and covariance matrix)
+
+       Int_t ndim = fNDim/3;
+       if(ndim>10) return 0.; // support only up to 10 dimensions
+
+       Int_t lambda = 1 + ndim + (ndim*(ndim+1)>>1);
+       Double_t fdfdp[66];
+       Int_t ipar = 0;
+       fdfdp[ipar++] = 1.;
+       for(int idim=0; idim<ndim; idim++){
+               fdfdp[ipar++] = point[idim];
+               for(int jdim=idim; jdim<ndim; jdim++) fdfdp[ipar++] = point[idim]*point[jdim];
+       }
+
+       // calculate estimation
+       result =0.; error = 0.;
+       for(int i=0; i<lambda; i++){
+               result += fdfdp[i]*(*fPar)(i);
+               for(int j=0; j<lambda; j++) error += fdfdp[i]*fdfdp[j]*(*fCov)(i,j);
+       }       
+       error = TMath::Sqrt(error);
+       
+       //printf("TKDNodeInfo::CookPDF() : %6.3f +- %6.3f\n", result, error);
+
+       return 0.;
+}
+
diff --git a/STAT/TKDNodeInfo.h b/STAT/TKDNodeInfo.h
new file mode 100644 (file)
index 0000000..9a522c9
--- /dev/null
@@ -0,0 +1,51 @@
+#ifndef ROOT_TKDNodeInfo
+#define ROOT_TKDNodeInfo
+
+#ifndef ROOT_TObject
+#include "TObject.h"
+#endif
+
+template <typename Value> class TVectorT;
+typedef struct TVectorT<Double_t> TVectorD;
+template <typename Value> class TMatrixT;
+typedef class TMatrixT<Double_t> TMatrixD;
+class TKDInterpolatorBase;
+class TKDNodeInfo : public TObject
+{
+friend class TKDInterpolatorBase;
+public:
+       TKDNodeInfo(const Int_t ndim = 0);
+       TKDNodeInfo(const TKDNodeInfo &);
+       TKDNodeInfo& operator=(const TKDNodeInfo &);
+       virtual  ~TKDNodeInfo();
+       Double_t  CookPDF(const Double_t *point, Double_t &result, Double_t &error);
+       inline Bool_t    Has(const Float_t *p);
+       void      Print();
+       void      Store(const TVectorD &par, const TMatrixD &cov);
+
+protected:
+       void      Build(const Int_t ndim);
+
+public:
+       Int_t     fNDim;          // 3 times data dimension
+       Float_t   *fData;         //[fNDim] node's data
+       Float_t   fVal[2];        // measured value for node
+       TMatrixD  *fCov;          // interpolator covariance matrix
+       TVectorD  *fPar;          // interpolator parameters
+
+       ClassDef(TKDNodeInfo, 1)  // node info for interpolator
+};
+
+//_____________________________________________________________________
+Bool_t TKDNodeInfo::Has(const Float_t *p)
+{
+       Int_t n = 0;
+       Int_t ndim = fNDim/3;
+       for(int id=0; id<ndim; id++) if(p[id]>=fData[ndim+2*id] && p[id]<fData[ndim+2*id+1]) n++;
+       if(n==ndim) return kTRUE;
+       return kFALSE;
+}
+
+
+#endif
+
diff --git a/STAT/TKDPDF.cxx b/STAT/TKDPDF.cxx
new file mode 100644 (file)
index 0000000..44f6b47
--- /dev/null
@@ -0,0 +1,184 @@
+#include "TKDPDF.h"
+#include "TKDNodeInfo.h"
+
+#include "TClonesArray.h"
+#include "TTree.h"
+#include "TH2.h"
+#include "TObjArray.h"
+#include "TObjString.h"
+#include "TBox.h"
+#include "TGraph.h"
+#include "TMarker.h"
+
+ClassImp(TKDPDF)
+
+
+
+//_________________________________________________________________
+TKDPDF::TKDPDF() :
+       TKDTreeIF()
+       ,TKDInterpolatorBase()
+{
+// Default constructor. To be used with care since in this case building
+// of data structure is completly left to the user responsability.
+}
+
+//_________________________________________________________________
+TKDPDF::TKDPDF(Int_t npoints, Int_t ndim, UInt_t bsize, Float_t **data) :
+       TKDTreeIF(npoints, ndim, bsize, data)
+       ,TKDInterpolatorBase(ndim)
+{
+// Wrapper constructor for the TKDTree.
+
+       Build();
+}
+
+
+//_________________________________________________________________
+TKDPDF::TKDPDF(TTree *t, const Char_t *var, const Char_t *cut, UInt_t bsize, Long64_t nentries, Long64_t firstentry) :
+       TKDTreeIF()
+       ,TKDInterpolatorBase()
+{
+// Alocate data from a tree. The variables which have to be analysed are
+// defined in the "var" parameter as a colon separated list. The format should
+// be identical to that used by TTree::Draw().
+//
+// 
+
+       TObjArray *vars = TString(var).Tokenize(":");
+       fNDim = vars->GetEntriesFast(); fNDimm = 2*fNDim;
+       fNSize = fNDim;
+       if(fNDim > 6/*kDimMax*/) Warning("TKDPDF(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", nentries, firstentry))){
+                       Warning("TKDPDF(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("TKDPDF(TTree*, const Char_t, const Char_t, UInt_t)", Form("Allocating %d data points in %d dimensions.", fNpoints, fNDim));
+                       fData = new Float_t*[fNDim];
+                       for(int idim=0; idim<fNDim; idim++) fData[idim] = new Float_t[fNpoints];
+                       kDataOwner = kTRUE;
+               }
+               v = t->GetV1();
+               for(int ip=0; ip<fNpoints; ip++) fData[idim][ip] = (Float_t)v[ip];
+       }
+       TKDTreeIF::Build();
+       Build();
+}
+
+//_________________________________________________________________
+TKDPDF::~TKDPDF()
+{
+}
+
+//_________________________________________________________________
+void TKDPDF::Build()
+{
+// Fill interpolator's data array i.e.
+//  - estimation points 
+//  - corresponding PDF values
+
+       fNTNodes = fNpoints/fBucketSize + ((fNpoints%fBucketSize)?1:0);/*TKDTreeIF::GetNTNodes();*/
+       if(!fBoundaries) MakeBoundaries();
+       fLambda = 1 + fNDim + (fNDim*(fNDim+1)>>1);
+       //printf("after MakeBoundaries() %d\n", memory());
+       
+       // allocate interpolation nodes
+       TKDInterpolatorBase::Build(fNTNodes);
+
+       TKDNodeInfo *node = 0x0;
+       Float_t *bounds = 0x0;
+       Int_t *indexPoints;
+       for(int inode=0, tnode = fNnodes; inode<fNTNodes-1; inode++, tnode++){
+               node = (TKDNodeInfo*)(*fTNodes)[inode];
+               node->fVal[0] =  Float_t(fBucketSize)/fNpoints;
+               bounds = GetBoundary(tnode);
+               for(int idim=0; idim<fNDim; idim++) node->fVal[0] /= (bounds[2*idim+1] - bounds[2*idim]);
+               node->fVal[1] =  node->fVal[0]/TMath::Sqrt(float(fBucketSize));
+               
+               indexPoints = GetPointsIndexes(tnode);
+               // loop points in this terminal node
+               for(int idim=0; idim<fNDim; idim++){
+                       node->fData[idim] = 0.;
+                       for(int ip = 0; ip<fBucketSize; ip++) node->fData[idim] += fData[idim][indexPoints[ip]];
+                       node->fData[idim] /= fBucketSize;
+               }
+               memcpy(&(node->fData[fNDim]), bounds, fNDimm*sizeof(Float_t));
+       }
+
+       // analyze last (incomplete) terminal node
+       Int_t counts = fNpoints%fBucketSize;
+       counts = counts ? counts : fBucketSize;
+       Int_t inode = fNTNodes - 1, tnode = inode + fNnodes;
+       node = (TKDNodeInfo*)(*fTNodes)[inode];
+       node->fVal[0] =  Float_t(counts)/fNpoints;
+       bounds = GetBoundary(tnode);
+       for(int idim=0; idim<fNDim; idim++) node->fVal[0] /= (bounds[2*idim+1] - bounds[2*idim]);
+       node->fVal[1] =  node->fVal[0]/TMath::Sqrt(float(counts));
+
+       // loop points in this terminal node
+       indexPoints = GetPointsIndexes(tnode);
+       for(int idim=0; idim<fNDim; idim++){
+               node->fData[idim] = 0.;
+               for(int ip = 0; ip<counts; ip++) node->fData[idim] += fData[idim][indexPoints[ip]];
+               node->fData[idim] /= counts;
+       }
+       memcpy(&(node->fData[fNDim]), bounds, fNDimm*sizeof(Float_t));
+
+       delete [] fBoundaries;
+       fBoundaries = 0x0;
+}
+
+
+//_________________________________________________________________
+void TKDPDF::DrawNode(Int_t tnode, UInt_t ax1, UInt_t ax2)
+{
+// Draw node "node" and the data points within.
+//
+// Observation:
+// This function creates some graphical objects
+// but don't delete it. Abusing this function may cause memory leaks !
+
+       if(tnode < 0 || tnode >= fNTNodes){
+               Warning("DrawNode()", Form("Terminal node %d outside defined range.", tnode));
+               return;
+       }
+
+       Int_t inode = tnode;
+       tnode += fNnodes;
+       // select zone of interest in the indexes array
+       Int_t *index = GetPointsIndexes(tnode);
+       Int_t nPoints = (tnode == 2*fNnodes) ? fNpoints%fBucketSize : fBucketSize;
+
+       // draw data points
+       TGraph *g = new TGraph(nPoints);
+       g->SetMarkerStyle(7);
+       for(int ip = 0; ip<nPoints; ip++) g->SetPoint(ip, fData[ax1][index[ip]], fData[ax2][index[ip]]);
+
+       // draw estimation point
+       TKDNodeInfo *node = (TKDNodeInfo*)(*fTNodes)[inode];
+       TMarker *m=new TMarker(node->fData[ax1], node->fData[ax2], 20);
+       m->SetMarkerColor(2);
+       m->SetMarkerSize(1.7);
+       
+       // draw node contour
+       Float_t *bounds = GetBoundary(tnode);
+       TBox *n = new TBox(bounds[2*ax1], bounds[2*ax2], bounds[2*ax1+1], bounds[2*ax2+1]);
+       n->SetFillStyle(0);
+
+       g->Draw("ap");
+       m->Draw();
+       n->Draw();
+       
+       return;
+}
+
diff --git a/STAT/TKDPDF.h b/STAT/TKDPDF.h
new file mode 100644 (file)
index 0000000..15f71a9
--- /dev/null
@@ -0,0 +1,59 @@
+#ifndef ROOT_TKDPDF
+#define ROOT_TKDPDF
+
+#ifndef ROOT_TKDTree
+#include "TKDTree.h"
+#endif
+
+#ifndef ROOT_TKDInterpolatorBase
+#include "TKDInterpolatorBase.h"
+#endif
+
+// Non parametric interpolation class based on local polinomial regression.
+// The class can be used to approximate PDF together with TKDTree or for
+// general regression when the data points are given.
+
+template <typename Index, typename Value> class TKDTree;
+typedef class TKDTree<Int_t, Float_t> TKDTreeIF;
+class TTree;
+class TLinearFitter;
+class TKDPDF : public TKDTreeIF, public TKDInterpolatorBase
+{
+public:
+       TKDPDF();
+       TKDPDF(TTree *t, const Char_t *var, const Char_t *cut = 0, UInt_t bsize = 100, Long64_t nentries = 1000000000, Long64_t firstentry = 0);
+       TKDPDF(Int_t npoints, Int_t ndim, UInt_t bsize, Float_t **data);
+       ~TKDPDF();
+
+       inline Bool_t     GetDataPoint(Int_t n, Float_t *p) const;
+       inline Int_t      GetNodeIndex(const Float_t *p);
+                                       void       DrawNode(Int_t tnode, UInt_t ax1=0, UInt_t ax2=1);
+
+private:
+       TKDPDF(const TKDPDF &);
+       TKDPDF& operator=(const TKDPDF &);
+               void       Build();
+
+                                       
+       ClassDef(TKDPDF, 1)   // data interpolator based on KD tree
+};
+
+
+//__________________________________________________________________
+Bool_t TKDPDF::GetDataPoint(Int_t n, Float_t *p) const
+{
+       if(n < 0 || n >= fNpoints) return kFALSE;
+       if(!fData) return kFALSE;
+               
+       for(int i=0; i<fNDim; i++) p[i] = fData[i][n];
+       return kTRUE;
+}
+
+//__________________________________________________________________
+Int_t  TKDPDF::GetNodeIndex(const Float_t *p)
+{
+       return FindNode(p) - fNnodes;
+}
+
+#endif
+
index 7f00d27..25315c7 100644 (file)
@@ -1,32 +1,18 @@
 #include "TKDSpline.h"
 
 
-
 ClassImp(TKDSpline)
 
-/////////////////////////////////////////////////////////////////////
-// Memory setup of protected data memebers
-// fRefPoints : evaluation point of PDF for each terminal node of underlying KD Tree.
-// | 1st terminal node (fNDim point coordinates) | 2nd terminal node (fNDim point coordinates) | ...
-//
-// 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) | ...
-/////////////////////////////////////////////////////////////////////
 
 //_________________________________________________________________
-TKDSpline::TKDSpline() : TKDInterpolator()
+TKDSpline::TKDSpline() :
+       TKDInterpolator()
 {
 }
 
 //_________________________________________________________________
-TKDSpline::TKDSpline(Int_t npoints, Int_t ndim, UInt_t bsize, Float_t **data) : TKDInterpolator(npoints, ndim, bsize, data)
-{
-       Build();
-}
-
-
-//_________________________________________________________________
-TKDSpline::TKDSpline(TTree *t, const Char_t *var, const Char_t *cut, UInt_t bsize) : TKDInterpolator(t, var, cut, bsize)
+TKDSpline::TKDSpline(Int_t npoints, Int_t ndim) :
+       TKDInterpolator(npoints, ndim)
 {
        Build();
 }
index c8c4ab8..a33741c 100644 (file)
@@ -9,8 +9,7 @@ class TKDSpline : public TKDInterpolator
 {
 public:
        TKDSpline();
-       TKDSpline(TTree *t, const Char_t *var, const Char_t *cut = 0, UInt_t bsize=100);
-       TKDSpline(Int_t npoints, Int_t ndim, UInt_t bsize, Float_t **data);
+       TKDSpline(Int_t npoints, Int_t ndim);
 
 private:
        void            Build();
index eebb080..775c602 100644 (file)
@@ -311,7 +311,7 @@ Bool_t      TKDTree<Index, Value>::FindNearestNeighbors(const Value *point, const Int
        // traverse tree
        UChar_t ax;
        Value   val, dif;
-       Int_t nAllNodes = fNnodes + GetNTNodes();
+       Int_t nAllNodes = fNnodes + fNpoints/fBucketSize + ((fNpoints%fBucketSize)?1:0);
        Int_t nodeStack[128], nodeIn[128];
        Index currentIndex = 0;
        nodeStack[0]   = inode;
@@ -745,7 +745,7 @@ void TKDTree<Index, Value>::MakeBoundaries(Value *range)
 
        if(range) memcpy(fRange, range, fNDimm*sizeof(Value));
        // total number of nodes including terminal nodes
-       Int_t totNodes = fNnodes + GetNTNodes();
+       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));
 
index 5ef6075..7d35d32 100644 (file)
@@ -28,7 +28,6 @@ public:
        inline UChar_t GetNodeAxis(Int_t id) const {return (id < 0 || id >= fNnodes) ? 0 : fAxis[id];}
        inline Value   GetNodeValue(Int_t id) const {return (id < 0 || id >= fNnodes) ? 0 : fValue[id];}
        inline Int_t   GetNNodes() const {return fNnodes;}
-       inline Int_t   GetNTNodes() 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);}
index ffc070c..e74c2d3 100644 (file)
@@ -1,5 +1,8 @@
 SRCS=  TKDTree.cxx \
-       TKDInterpolator.cxx \
+  TKDInterpolatorBase.cxx \
+  TKDNodeInfo.cxx \
+  TKDPDF.cxx \
+  TKDInterpolator.cxx \
   TKDSpline.cxx \
   AliTMinuitToolkit.cxx \
   TStatToolkit.cxx