From: abercuci Date: Thu, 11 Dec 2008 19:42:13 +0000 (+0000) Subject: TKDTree class now in ROOT X-Git-Url: http://git.uio.no/git/?a=commitdiff_plain;h=50c4eb3a098aa07518004665419cb2cbc9609ad8;p=u%2Fmrichter%2FAliRoot.git TKDTree class now in ROOT --- diff --git a/STAT/CMake_libSTAT.txt b/STAT/CMake_libSTAT.txt index 91e7bc2d32b..fa7eaade5da 100644 --- a/STAT/CMake_libSTAT.txt +++ b/STAT/CMake_libSTAT.txt @@ -1,7 +1,6 @@ # -*- mode: cmake -*- set(SRCS - TKDTree.cxx TKDInterpolatorBase.cxx TKDNodeInfo.cxx TKDPDF.cxx diff --git a/STAT/STATLinkDef.h b/STAT/STATLinkDef.h index 3fa52539bcb..5b09bacb7a6 100644 --- a/STAT/STATLinkDef.h +++ b/STAT/STATLinkDef.h @@ -7,12 +7,6 @@ #pragma link off all classes; #pragma link off all functions; -#pragma link C++ class TKDTree+; -#pragma link C++ class TKDTree+; -#pragma link C++ typedef TKDTreeID; -#pragma link C++ typedef TKDTreeIF; - - #pragma link C++ class TKDInterpolatorBase+; #pragma link C++ class TKDNodeInfo+; #pragma link C++ class TKDPDF+; diff --git a/STAT/TKDInterpolatorBase.cxx b/STAT/TKDInterpolatorBase.cxx index 1999b959c10..3ae38700b0a 100644 --- a/STAT/TKDInterpolatorBase.cxx +++ b/STAT/TKDInterpolatorBase.cxx @@ -31,17 +31,17 @@ ClassImp(TKDInterpolatorBase) //_________________________________________________________________ TKDInterpolatorBase::TKDInterpolatorBase(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) + 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. @@ -50,46 +50,46 @@ TKDInterpolatorBase::TKDInterpolatorBase(Int_t dim) : //_________________________________________________________________ void TKDInterpolatorBase::Build(Int_t n) { - // allocate memory for data + // allocate memory for data - if(fTNodes) delete fTNodes; - fNTNodes = n; - fTNodes = new TClonesArray("TKDNodeInfo", fNTNodes); - for(int in=0; in fNTNodes) return kFALSE; + if(inode < 0 || inode > fNTNodes) return kFALSE; - TKDNodeInfo *node = (TKDNodeInfo*)(*fTNodes)[inode]; - coord = &(node->Data()[0]); - val = node->Val()[0]; - err = node->Val()[1]; - return kTRUE; + TKDNodeInfo *node = (TKDNodeInfo*)(*fTNodes)[inode]; + coord = &(node->Data()[0]); + val = node->Val()[0]; + err = node->Val()[1]; + return kTRUE; } //_________________________________________________________________ TKDNodeInfo* TKDInterpolatorBase::GetNodeInfo(Int_t inode) const { - if(!fTNodes || inode >= fNTNodes) return 0x0; - return (TKDNodeInfo*)(*fTNodes)[inode]; + if(!fTNodes || inode >= fNTNodes) return 0x0; + return (TKDNodeInfo*)(*fTNodes)[inode]; } @@ -98,17 +98,17 @@ 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; iPrint(); - } + 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; iPrint(); + } } //_________________________________________________________________ @@ -120,119 +120,116 @@ Double_t TKDInterpolatorBase::Eval(const Double_t *point, Double_t &result, Doub // // 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; idimCov() && !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; idData()[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; idimFindNearestNeighbors(pointF, npoints+1, index, dist)){ - Error("Eval()", Form("Failed retriving %d neighbours for point:", npoints)); - for(int idim=0; idimClearPoints(); - TKDNodeInfo *tnode = 0x0; - for(int in=0; inPrint(); - if(fStatus&1){ // INT - Float_t *bounds = &(tnode->Data()[fNSize]); - ipar = 0; - for(int idim=0; idimData()[0]); - ipar = 0; - for(int idim=0; idimCov() && !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; idData()[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 ipar, // local looping variable + npoints = Int_t((1.+fAlpha)*fLambda); // number of data points used for interpolation + Int_t *index = new Int_t[2*npoints]; // indexes of NN + Float_t *dist = new Float_t[2*npoints], // distances of NN + d, // NN normalized distance + w0, // work + w; // tri-cubic weight function + + do{ + // find nearest neighbors + for(int idim=0; idimFindNearestNeighbors(pointF, npoints+1, index, dist); + // add points to fitter + fFitter->ClearPoints(); + TKDNodeInfo *tnode = 0x0; + for(int in=0; inPrint(); + if(fStatus&1){ // INT + Float_t *bounds = &(tnode->Data()[fNSize]); + ipar = 0; + for(int idim=0; idimData()[0]); + ipar = 0; + for(int idim=0; idimVal()[0], tnode->Val()[1]/w, tnode->Val()[1], w); - fFitter->AddPoint(fBuffer, tnode->Val()[0], tnode->Val()[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; idimAddPoint(fBuffer, tnode->Val()[0], tnode->Val()[1]/w); + } + npoints += 4; + } while(fFitter->Eval()); + delete [] index; + delete [] dist; + + // 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; idimGetXaxis()->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])->Data()[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->Data()[ax1], node->Data()[ax2]); - } - ref->Draw("p"); - return; + + 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])->Data()[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->Data()[ax1], node->Data()[ax2]); + } + ref->Draw("p"); + return; } //__________________________________________________________________ void TKDInterpolatorBase::SetInterpolationMethod(Bool_t on) { // Set interpolation bit to "on". - - if(on) fStatus += fStatus&1 ? 0 : 1; - else fStatus += fStatus&1 ? -1 : 0; + + if(on) fStatus += fStatus&1 ? 0 : 1; + else fStatus += fStatus&1 ? -1 : 0; } @@ -292,16 +289,16 @@ void TKDInterpolatorBase::SetInterpolationMethod(Bool_t on) void TKDInterpolatorBase::SetStore(Bool_t on) { // Set store bit to "on" - - if(on) fStatus += fStatus&2 ? 0 : 2; - else fStatus += fStatus&2 ? -2 : 0; + + if(on) fStatus += fStatus&2 ? 0 : 2; + else fStatus += fStatus&2 ? -2 : 0; } //_________________________________________________________________ void TKDInterpolatorBase::SetWeights(Bool_t on) { // Set weights bit to "on" - - if(on) fStatus += fStatus&4 ? 0 : 4; - else fStatus += fStatus&4 ? -4 : 0; + + if(on) fStatus += fStatus&4 ? 0 : 4; + else fStatus += fStatus&4 ? -4 : 0; } diff --git a/STAT/TKDInterpolatorBase.h b/STAT/TKDInterpolatorBase.h index ae7097b7bb3..2674cf2c72f 100644 --- a/STAT/TKDInterpolatorBase.h +++ b/STAT/TKDInterpolatorBase.h @@ -26,49 +26,49 @@ class TKDNodeInfo; class TKDInterpolatorBase { public: - TKDInterpolatorBase(Int_t size = 0); - virtual ~TKDInterpolatorBase(); + TKDInterpolatorBase(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(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;} + 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(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(Float_t a){if(a>0.) fAlpha = a;} - void SetInterpolationMethod(Bool_t on = kTRUE); - void SetStore(Bool_t on = kTRUE); - void SetWeights(Bool_t on = kTRUE); + 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(Float_t a){if(a>0.) fAlpha = a;} + void SetInterpolationMethod(Bool_t on = kTRUE); + void SetStore(Bool_t on = kTRUE); + void SetWeights(Bool_t on = kTRUE); protected: - virtual void Build(Int_t nnodes); + virtual void Build(Int_t nnodes); private: - TKDInterpolatorBase(const TKDInterpolatorBase &); - TKDInterpolatorBase& operator=(const TKDInterpolatorBase &); + 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 + 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 *fKDhelper; //! kNN finder - TLinearFitter *fFitter; //! linear fitter + 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 *fKDhelper; //! kNN finder + TLinearFitter *fFitter; //! linear fitter - ClassDef(TKDInterpolatorBase, 1) // data interpolator based on KD tree + ClassDef(TKDInterpolatorBase, 1) // data interpolator based on KD tree }; diff --git a/STAT/TKDPDF.cxx b/STAT/TKDPDF.cxx index c8fa07e6460..72bf23ad8e7 100644 --- a/STAT/TKDPDF.cxx +++ b/STAT/TKDPDF.cxx @@ -16,8 +16,8 @@ ClassImp(TKDPDF) //_________________________________________________________________ TKDPDF::TKDPDF() : - TKDTreeIF() - ,TKDInterpolatorBase() + TKDTreeIF() + ,TKDInterpolatorBase() { // Default constructor. To be used with care since in this case building // of data structure is completly left to the user responsability. @@ -25,19 +25,19 @@ TKDPDF::TKDPDF() : //_________________________________________________________________ TKDPDF::TKDPDF(Int_t npoints, Int_t ndim, UInt_t bsize, Float_t **data) : - TKDTreeIF(npoints, ndim, bsize, data) - ,TKDInterpolatorBase(ndim) + TKDTreeIF(npoints, ndim, bsize, data) + ,TKDInterpolatorBase(ndim) { // Wrapper constructor for the TKDTree. - Build(); + Build(); } //_________________________________________________________________ TKDPDF::TKDPDF(TTree *t, const Char_t *var, const Char_t *cut, UInt_t bsize, Long64_t nentries, Long64_t firstentry) : - TKDTreeIF() - ,TKDInterpolatorBase() + 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 @@ -45,34 +45,34 @@ TKDPDF::TKDPDF(TTree *t, const Char_t *var, const Char_t *cut, UInt_t bsize, Lon // // - 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; idimDraw(((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; idimGetV1(); - for(int ip=0; ipGetEntriesFast(); 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; idimDraw(((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; idimGetV1(); + for(int ip=0; ip>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; inodeVal()[0] = Float_t(fBucketSize)/fNpoints; - bounds = GetBoundary(tnode); - for(int idim=0; idimVal()[0] /= (bounds[2*idim+1] - bounds[2*idim]); - node->Val()[1] = node->Val()[0]/TMath::Sqrt(float(fBucketSize)); - - indexPoints = GetPointsIndexes(tnode); - // loop points in this terminal node - for(int idim=0; idimData()[idim] = 0.; - for(int ip = 0; ipData()[idim] += fData[idim][indexPoints[ip]]; - node->Data()[idim] /= fBucketSize; - } - memcpy(&(node->Data()[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->Val()[0] = Float_t(counts)/fNpoints; - bounds = GetBoundary(tnode); - for(int idim=0; idimVal()[0] /= (bounds[2*idim+1] - bounds[2*idim]); - node->Val()[1] = node->Val()[0]/TMath::Sqrt(float(counts)); - - // loop points in this terminal node - indexPoints = GetPointsIndexes(tnode); - for(int idim=0; idimData()[idim] = 0.; - for(int ip = 0; ipData()[idim] += fData[idim][indexPoints[ip]]; - node->Data()[idim] /= counts; - } - memcpy(&(node->Data()[fNDim]), bounds, fNDimm*sizeof(Float_t)); - - delete [] fBoundaries; - fBoundaries = 0x0; + 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; inodeVal()[0] = Float_t(fBucketSize)/fNPoints; + bounds = GetBoundary(tnode); + for(int idim=0; idimVal()[0] /= (bounds[2*idim+1] - bounds[2*idim]); + node->Val()[1] = node->Val()[0]/TMath::Sqrt(float(fBucketSize)); + + indexPoints = GetPointsIndexes(tnode); + // loop points in this terminal node + for(int idim=0; idimData()[idim] = 0.; + for(int ip = 0; ipData()[idim] += fData[idim][indexPoints[ip]]; + node->Data()[idim] /= fBucketSize; + } + memcpy(&(node->Data()[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->Val()[0] = Float_t(counts)/fNPoints; + bounds = GetBoundary(tnode); + for(int idim=0; idimVal()[0] /= (bounds[2*idim+1] - bounds[2*idim]); + node->Val()[1] = node->Val()[0]/TMath::Sqrt(float(counts)); + + // loop points in this terminal node + indexPoints = GetPointsIndexes(tnode); + for(int idim=0; idimData()[idim] = 0.; + for(int ip = 0; ipData()[idim] += fData[idim][indexPoints[ip]]; + node->Data()[idim] /= counts; + } + memcpy(&(node->Data()[fNDim]), bounds, fNDimm*sizeof(Float_t)); + + delete [] fBoundaries; + fBoundaries = 0x0; } @@ -148,37 +148,37 @@ void TKDPDF::DrawNode(Int_t tnode, UInt_t ax1, UInt_t ax2) // 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; ipSetPoint(ip, fData[ax1][index[ip]], fData[ax2][index[ip]]); - - // draw estimation point - TKDNodeInfo *node = (TKDNodeInfo*)(*fTNodes)[inode]; - TMarker *m=new TMarker(node->Data()[ax1], node->Data()[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; + 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; ipSetPoint(ip, fData[ax1][index[ip]], fData[ax2][index[ip]]); + + // draw estimation point + TKDNodeInfo *node = (TKDNodeInfo*)(*fTNodes)[inode]; + TMarker *m=new TMarker(node->Data()[ax1], node->Data()[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 index 846e3e05340..8ec6a384ce0 100644 --- a/STAT/TKDPDF.h +++ b/STAT/TKDPDF.h @@ -1,14 +1,14 @@ #ifndef ROOT_TKDPDF #define ROOT_TKDPDF -#ifndef ROOT_TKDTree -#include "TKDTree.h" -#endif - #ifndef ROOT_TKDInterpolatorBase #include "TKDInterpolatorBase.h" #endif +#ifndef ROOT_TKDTree +#include "TKDTree.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. @@ -42,7 +42,7 @@ private: //__________________________________________________________________ Bool_t TKDPDF::GetDataPoint(Int_t n, Float_t *p) const { - if(n < 0 || n >= fNpoints) return kFALSE; + if(n < 0 || n >= fNPoints) return kFALSE; if(!fData) return kFALSE; for(int i=0; i