STAT/TKDInterpolatorBase.cxx

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