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