STAT/TKDInterpolatorBase.cxx

   1 #include "TKDInterpolatorBase.h"
   2 #include "TKDNodeInfo.h"
   3 #include "TKDTree.h"
   4
   5 #include "TROOT.h"
   6 #include "TClonesArray.h"
   7 #include "TLinearFitter.h"
   8 #include "TTree.h"
   9 #include "TH2.h"
  10 #include "TObjArray.h"
  11 #include "TObjString.h"
  12 #include "TMath.h"
  13 #include "TVectorD.h"
  14 #include "TMatrixD.h"
  15
  16 ClassImp(TKDInterpolatorBase)
  17
  18 /////////////////////////////////////////////////////////////////////
  19 // Memory setup of protected data members
  20 // fRefPoints : evaluation point of PDF for each terminal node of underlying KD Tree.
  21 // | 1st terminal node (fNDim point coordinates) | 2nd terminal node (fNDim point coordinates) | ...
  22 //
  23 // fRefValues : evaluation value/error of PDF for each terminal node of underlying KD Tree.
  24 // | 1st terminal node (value) | 2nd terminal node (value) | ... | 1st terminal node (error) | 2nd terminal node (error) | ...
  25 //
  26 // status = |0|0|0|0|0|1(tri-cubic weights)|1(STORE)|1 INT(0 COG )|
  27 /////////////////////////////////////////////////////////////////////
  28
  29
  30 //_________________________________________________________________
  31 TKDInterpolatorBase::TKDInterpolatorBase(Int_t dim) :
  32   fNSize(dim)
  33   ,fNTNodes(0)
  34   ,fTNodes(0x0)
  35   ,fTNodesDraw(0x0)
  36   ,fStatus(0)
  37   ,fLambda(1 + dim + (dim*(dim+1)>>1))
  38   ,fDepth(-1)
  39   ,fAlpha(.5)
  40   ,fRefPoints(0x0)
  41   ,fBuffer(0x0)
  42   ,fKDhelper(0x0)
  43   ,fFitter(0x0)
  44 {
  45 // Default constructor. To be used with care since in this case building
  46 // of data structure is completly left to the user responsability.
  47   UseWeights();
  48 }
  49
  50 //_________________________________________________________________
  51 void    TKDInterpolatorBase::Build(Int_t n)
  52 {
  53   // allocate memory for data
  54
  55   if(fTNodes) delete fTNodes;
  56   fNTNodes = n;
  57   // check granularity
  58   if(Int_t((1.+fAlpha)*fLambda) > fNTNodes){
  59     Warning("TKDInterpolatorBase::Build()", Form("Minimum number of points [%d] needed for interpolation exceeds number of evaluation points [%d]. Please increase granularity.", Int_t((1.+fAlpha)*fLambda), fNTNodes));
  60   }
  61   fTNodes = new TClonesArray("TKDNodeInfo", fNTNodes);
  62   for(int in=0; in<fNTNodes; in++) new ((*fTNodes)[in]) TKDNodeInfo(fNSize);
  63 }
  64
  65 //_________________________________________________________________
  66 TKDInterpolatorBase::~TKDInterpolatorBase()
  67 {
  68   if(fFitter) delete fFitter;
  69   if(fKDhelper) delete fKDhelper;
  70   if(fBuffer) delete [] fBuffer;
  71
  72   if(fRefPoints){
  73     for(int idim=0; idim<fNSize; idim++) delete [] fRefPoints[idim] ;
  74     delete [] fRefPoints;
  75   }
  76   if(fTNodes){
  77     fTNodes->Delete();
  78     delete fTNodes;
  79   }
  80   delete [] fTNodesDraw;
  81
  82   TH2 *h2=0x0;
  83   if((h2 = (TH2S*)gROOT->FindObject("hKDnodes"))) delete h2;
  84 }
  85
  86
  87 //__________________________________________________________________
  88 Bool_t  TKDInterpolatorBase::GetCOGPoint(Int_t inode, Float_t *&coord, Float_t &val, Float_t &err) const
  89 {
  90   if(inode < 0 || inode > fNTNodes) return kFALSE;
  91
  92   TKDNodeInfo *node = (TKDNodeInfo*)(*fTNodes)[inode];
  93   coord = &(node->Data()[0]);
  94   val = node->Val()[0];
  95   err = node->Val()[1];
  96   return kTRUE;
  97 }
  98
  99 //_________________________________________________________________
 100 TKDNodeInfo* TKDInterpolatorBase::GetNodeInfo(Int_t inode) const
 101 {
 102   if(!fTNodes || inode >= fNTNodes) return 0x0;
 103   return (TKDNodeInfo*)(*fTNodes)[inode];
 104 }
 105
 106 //_________________________________________________________________
 107 Bool_t TKDInterpolatorBase::GetRange(Int_t ax, Float_t &min, Float_t &max) const
 108 {
 109   if(!fTNodes) return kFALSE;
 110   Int_t ndim = ((TKDNodeInfo*)(*fTNodes)[0])->GetDimension();
 111   if(ax<0 || ax>=ndim){
 112     min=0.; max=0.;
 113     return kFALSE;
 114   }
 115   min=1.e10; max=-1.e10;
 116   Float_t axmin, axmax;
 117   for(Int_t in=fNTNodes; in--; ){
 118     TKDNodeInfo *node = (TKDNodeInfo*)((*fTNodes)[in]);
 119     node->GetBoundary(ax, axmin, axmax);
 120     if(axmin<min) min = axmin;
 121     if(axmax>max) max = axmax;
 122   }
 123
 124   return kTRUE;
 125 }
 126
 127 //__________________________________________________________________
 128 void TKDInterpolatorBase::GetStatus(Option_t *opt)
 129 {
 130 // Prints the status of the interpolator
 131
 132   printf("Interpolator Status[%d] :\n", fStatus);
 133   printf("  Dim    : %d [%d]\n", fNSize, fLambda);
 134   printf("  Method : %s\n", UseCOG() ? "COG" : "INT");
 135   printf("  Store  : %s\n", HasStore() ? "YES" : "NO");
 136   printf("  Weights: %s\n", UseWeights() ? "YES" : "NO");
 137
 138   if(strcmp(opt, "all") != 0 ) return;
 139   printf("fNTNodes %d\n", fNTNodes);        //Number of evaluation data points
 140   for(int i=0; i<fNTNodes; i++){
 141     TKDNodeInfo *node = (TKDNodeInfo*)(*fTNodes)[i];
 142     printf("%d ", i); node->Print();
 143   }
 144 }
 145
 146 //_________________________________________________________________
 147 Double_t TKDInterpolatorBase::Eval(const Double_t *point, Double_t &result, Double_t &error, Bool_t force)
 148 {
 149 // Evaluate PDF for "point". The result is returned in "result" and error in "error". The function returns the chi2 of the fit.
 150 //
 151 // Observations:
 152 //
 153 // 1. The default method used for interpolation is kCOG.
 154 // 2. The initial number of neighbors used for the estimation is set to Int(alpha*fLambda) (alpha = 1.5)
 155
 156   Float_t pointF[50]; // local Float_t conversion for "point"
 157   for(int idim=0; idim<fNSize; idim++) pointF[idim] = (Float_t)point[idim];
 158   Int_t nodeIndex = GetNodeIndex(pointF);
 159   if(nodeIndex<0){
 160     Error("TKDInterpolatorBase::Eval()", "Can not retrive node for data point.");
 161     result = 0.;
 162     error = 1.E10;
 163     return 0.;
 164   }
 165   TKDNodeInfo *node = (TKDNodeInfo*)(*fTNodes)[nodeIndex];
 166   if(node->Cov() && !force) return node->CookPDF(point, result, error);
 167
 168   // Allocate memory
 169   if(!fBuffer) fBuffer = new Double_t[2*fLambda];
 170   if(!fKDhelper){
 171     fRefPoints = new Float_t*[fNSize];
 172     for(int id=0; id<fNSize; id++){
 173       fRefPoints[id] = new Float_t[fNTNodes];
 174       for(int in=0; in<fNTNodes; in++) fRefPoints[id][in] = ((TKDNodeInfo*)(*fTNodes)[in])->Data()[id];
 175     }
 176     Info("TKDInterpolatorBase::Eval()", Form("Build TKDTree(%d, %d, %d)", fNTNodes, fNSize, kNhelper));
 177     fKDhelper = new TKDTreeIF(fNTNodes, fNSize, kNhelper, fRefPoints);
 178     fKDhelper->Build();
 179     fKDhelper->MakeBoundariesExact();
 180   }
 181   if(!fFitter) fFitter = new TLinearFitter(fLambda, Form("hyp%d", fLambda-1));
 182
 183   // generate parabolic for nD
 184   //Float_t alpha = Float_t(2*lambda + 1) / fNTNodes; // the bandwidth or smoothing parameter
 185   //Int_t npoints = Int_t(alpha * fNTNodes);
 186   //printf("Params : %d NPoints %d\n", lambda, npoints);
 187   // prepare workers
 188
 189   Int_t ipar,    // local looping variable
 190         npoints_new = Int_t((1.+fAlpha)*fLambda),
 191         npoints(0); // number of data points used for interpolation
 192   Int_t *index = new Int_t[2*npoints_new];  // indexes of NN
 193   Float_t *dist = new Float_t[2*npoints_new], // distances of NN
 194           d,     // NN normalized distance
 195           w0,    // work
 196           w;     // tri-cubic weight function
 197
 198   Bool_t kDOWN = kFALSE;
 199   do{
 200     if(npoints){
 201       Info("TKDInterpolatorBase::Eval()", Form("Interpolation failed. Trying to increase the number of interpolation points from %d to %d.", npoints, npoints_new));
 202     }
 203     if(npoints == npoints_new){
 204       Error("TKDInterpolatorBase::Eval()", Form("Interpolation failed and number of interpolation points (%d) Can not be increased further.", npoints));
 205       result = 0.;
 206       error = 1.E10;
 207       return 0.;
 208     } else npoints = npoints_new;
 209     if(npoints > fNTNodes){
 210       Warning("TKDInterpolatorBase::Eval()", Form("The number of interpolation points requested (%d) exceeds number of PDF values (%d). Downscale.", npoints, fNTNodes));
 211       npoints = fNTNodes;
 212       kDOWN = kTRUE;
 213     }
 214
 215     // find nearest neighbors
 216     for(int idim=0; idim<fNSize; idim++) pointF[idim] = (Float_t)point[idim];
 217     fKDhelper->FindNearestNeighbors(pointF, npoints+1, index, dist);
 218
 219     // add points to fitter
 220     fFitter->ClearPoints();
 221     TKDNodeInfo *tnode = 0x0;
 222     for(int in=0; in<npoints; in++){
 223       tnode = (TKDNodeInfo*)(*fTNodes)[index[in]];
 224       //tnode->Print();
 225       if(UseCOG()){ // COG
 226         Float_t *p = &(tnode->Data()[0]);
 227         ipar = 0;
 228         for(int idim=0; idim<fNSize; idim++){
 229           fBuffer[ipar++] = p[idim];
 230           for(int jdim=idim; jdim<fNSize; jdim++) fBuffer[ipar++] = p[idim]*p[jdim];
 231         }
 232       } else { // INT
 233         Float_t *bounds = &(tnode->Data()[fNSize]);
 234         ipar = 0;
 235         for(int idim=0; idim<fNSize; idim++){
 236           fBuffer[ipar++] = .5*(bounds[2*idim] + bounds[2*idim+1]);
 237           fBuffer[ipar++] = (bounds[2*idim]*bounds[2*idim] + bounds[2*idim] * bounds[2*idim+1] + bounds[2*idim+1] * bounds[2*idim+1])/3.;
 238           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;
 239         }
 240       }
 241
 242       // calculate tri-cubic weighting function
 243       if(UseWeights()){
 244         d = dist[in]/dist[npoints];
 245         w0 = (1. - d*d*d); w = w0*w0*w0;
 246         if(w<1.e-30) continue;
 247       } else w = 1.;
 248
 249 //       printf("%2d d[%f] w[%f] x[", index[in], d, w);
 250 //       for(int idim=0; idim<fLambda-1; idim++) printf("%f ", fBuffer[idim]);
 251 //       printf("]\n");  printf("v[%f +- %f] (%f, %f)\n", tnode->Val()[0], tnode->Val()[1]/w, tnode->Val()[1], w);
 252       fFitter->AddPoint(fBuffer, tnode->Val()[0], tnode->Val()[1]/w);
 253     }
 254     npoints_new = npoints+ (kDOWN ? 0 : kdN);
 255   } while(fFitter->Eval());
 256   delete [] index;
 257   delete [] dist;
 258
 259   // retrive fitter results
 260   TMatrixD cov(fLambda, fLambda);
 261   TVectorD par(fLambda);
 262   fFitter->GetCovarianceMatrix(cov);
 263   fFitter->GetParameters(par);
 264   Double_t chi2 = fFitter->GetChisquare()/(npoints - 4 - fLambda);
 265
 266   // store results
 267   if(HasStore()) node->Store(par, cov);
 268
 269   // Build df/dpi|x values
 270   Double_t *fdfdp = &fBuffer[fLambda];
 271   ipar = 0;
 272   fdfdp[ipar++] = 1.;
 273   for(int idim=0; idim<fNSize; idim++){
 274     fdfdp[ipar++] = point[idim];
 275     for(int jdim=idim; jdim<fNSize; jdim++) fdfdp[ipar++] = point[idim]*point[jdim];
 276   }
 277
 278   // calculate estimation
 279   result =0.; error = 0.;
 280   for(int i=0; i<fLambda; i++){
 281     result += fdfdp[i]*par(i);
 282     for(int j=0; j<fLambda; j++) error += fdfdp[i]*fdfdp[j]*cov(i,j);
 283   }
 284   error = TMath::Sqrt(error);
 285
 286   return chi2;
 287 }
 288
 289 //_________________________________________________________________
 290 void TKDInterpolatorBase::DrawProjection(UInt_t ax1, UInt_t ax2)
 291 {
 292 // Draw nodes structure projected on plane "ax1:ax2". The parameter
 293 // "depth" specifies the bucket size per node. If depth == -1 draw only
 294 // terminal nodes and evaluation points (default -1 i.e. bucket size per node equal bucket size specified by the user)
 295 //
 296
 297   Float_t ax1min, ax1max, ax2min, ax2max;
 298   GetRange(ax1, ax1min, ax1max);
 299   GetRange(ax2, ax2min, ax2max);
 300   TH2 *h2 = 0x0;
 301   if(!(h2 = (TH2S*)gROOT->FindObject("hKDnodes"))){
 302     h2 = new TH2S("hKDnodes", "", 100, ax1min, ax1max, 100, ax2min, ax2max);
 303   }
 304   h2->GetXaxis()->SetRangeUser(ax1min, ax1max);
 305   h2->GetXaxis()->SetTitle(Form("x_{%d}", ax1));
 306   h2->GetYaxis()->SetRangeUser(ax2min, ax2max);
 307   h2->GetYaxis()->SetTitle(Form("x_{%d}", ax2));
 308   h2->Draw();
 309
 310
 311   if(!fTNodesDraw) fTNodesDraw = new TKDNodeInfo::TKDNodeDraw[fNTNodes];
 312   TKDNodeInfo::TKDNodeDraw *box = 0x0;
 313   for(Int_t in=fNTNodes; in--; ){
 314     box = &(fTNodesDraw[in]);
 315     box->SetNode((TKDNodeInfo*)((*fTNodes)[in]), fNSize, ax1, ax2);
 316     box->Draw();
 317   }
 318
 319   return;
 320 }
 321
 322 //_________________________________________________________________
 323 void TKDInterpolatorBase::SetAlpha(Float_t a)
 324 {
 325   if(a<0.5){
 326     Warning("TKDInterpolatorBase::SetAlpha()", "The scale parameter has to be larger than 0.5");
 327     fAlpha = 0.5;
 328     return;
 329   }
 330   // check value
 331   if(Int_t((a+1.)*fLambda) > fNTNodes){
 332     fAlpha = TMath::Max(0.5, Float_t(fNTNodes)/fLambda-1.);
 333     Warning("TKDInterpolatorBase::SetAlpha()", Form("Interpolation neighborhood  exceeds number of evaluation points. Downscale alpha to %f", fAlpha));
 334     printf("n[%d] nodes[%d]\n", Int_t((fAlpha+1.)*fLambda), fNTNodes);
 335     return;
 336   }
 337   fAlpha = a;
 338   return;
 339 }
 340