From a3408ed3d4c63c1758fe76fa0e6d0d5b8c268038 Mon Sep 17 00:00:00 2001 From: abercuci Date: Thu, 25 Oct 2007 14:04:29 +0000 Subject: [PATCH] rename Interpolator to PDF, add new class TKDInterpolator and modify inheritance tree --- STAT/StatLinkDef.h | 4 +- STAT/TKDInterpolator.cxx | 529 +++-------------------------------- STAT/TKDInterpolator.h | 104 +------ STAT/TKDInterpolatorBase.cxx | 307 ++++++++++++++++++++ STAT/TKDInterpolatorBase.h | 76 +++++ STAT/TKDNodeInfo.cxx | 152 ++++++++++ STAT/TKDNodeInfo.h | 51 ++++ STAT/TKDPDF.cxx | 184 ++++++++++++ STAT/TKDPDF.h | 59 ++++ STAT/TKDSpline.cxx | 22 +- STAT/TKDSpline.h | 3 +- STAT/TKDTree.cxx | 4 +- STAT/TKDTree.h | 1 - STAT/libSTAT.pkg | 5 +- 14 files changed, 891 insertions(+), 610 deletions(-) create mode 100644 STAT/TKDInterpolatorBase.cxx create mode 100644 STAT/TKDInterpolatorBase.h create mode 100644 STAT/TKDNodeInfo.cxx create mode 100644 STAT/TKDNodeInfo.h create mode 100644 STAT/TKDPDF.cxx create mode 100644 STAT/TKDPDF.h diff --git a/STAT/StatLinkDef.h b/STAT/StatLinkDef.h index dcbb2583831..3fa52539bcb 100644 --- a/STAT/StatLinkDef.h +++ b/STAT/StatLinkDef.h @@ -13,8 +13,10 @@ #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+; diff --git a/STAT/TKDInterpolator.cxx b/STAT/TKDInterpolator.cxx index 1d643368b7d..b9d0c1fd509 100644 --- a/STAT/TKDInterpolator.cxx +++ b/STAT/TKDInterpolator.cxx @@ -1,535 +1,82 @@ #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; idimGetEntriesFast(); 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; idimDraw(((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; idimGetV1(); - for(int ip=0; ip>1); - //printf("after MakeBoundaries() %d\n", memory()); - - // allocate memory for data - fTNodes = new TKDNodeInfo[fNTNodes]; - for(int in=0; inMakeBoundaries(); - } - 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; idimFindNearestNeighbors(pointF, npoints+1, index, dist)){ - Error("Eval()", Form("Failed retriving %d neighbours for point:", npoints)); - for(int idim=0; idimClearPoints(); - TKDNodeInfo *node = 0x0; - for(int in=0; infRefPoint)); - ipar = 0; - for(int idim=0; idimfRefPoint; - ipar = 0; - for(int idim=0; idimAddPoint(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> 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; iPrint(); + 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; ipSetPoint(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; -} diff --git a/STAT/TKDInterpolator.h b/STAT/TKDInterpolator.h index b8913bf18d9..94b47837c97 100644 --- a/STAT/TKDInterpolator.h +++ b/STAT/TKDInterpolator.h @@ -1,114 +1,30 @@ #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 class TVectorT; -typedef struct TVectorT TVectorD; -template class TMatrixT; -typedef class TMatrixT TMatrixD; -template class TKDTree; -typedef class TKDTree 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= fNpoints) return kFALSE; - if(!fData) return kFALSE; - - for(int i=0; i>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) 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; iPrint(); + } +} + +//_________________________________________________________________ +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; idimfCov && !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; idfData[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->fData[fNSize]); + ipar = 0; + for(int idim=0; idimfData[0]); + ipar = 0; + for(int idim=0; idimfVal[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; 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])->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 index 00000000000..5b591703655 --- /dev/null +++ b/STAT/TKDInterpolatorBase.h @@ -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 class TVectorT; +typedef struct TVectorT TVectorD; +template class TMatrixT; +typedef class TMatrixT TMatrixD; +template class TKDTree; +typedef class TKDTree 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 *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 index 00000000000..a4069520624 --- /dev/null +++ b/STAT/TKDNodeInfo.cxx @@ -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; idim10) 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 class TVectorT; +typedef struct TVectorT TVectorD; +template class TMatrixT; +typedef class TMatrixT 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=fData[ndim+2*id] && p[id]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; 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; inodefVal[0] = Float_t(fBucketSize)/fNpoints; + bounds = GetBoundary(tnode); + for(int idim=0; idimfVal[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; idimfData[idim] = 0.; + for(int ip = 0; ipfData[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; idimfVal[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; idimfData[idim] = 0.; + for(int ip = 0; ipfData[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; ipSetPoint(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 index 00000000000..15f71a9971a --- /dev/null +++ b/STAT/TKDPDF.h @@ -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 class TKDTree; +typedef class TKDTree 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::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::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)); diff --git a/STAT/TKDTree.h b/STAT/TKDTree.h index 5ef60752e9d..7d35d3299bb 100644 --- a/STAT/TKDTree.h +++ b/STAT/TKDTree.h @@ -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<