From Philippe & Laurent: new variant of MUON visualization.
[u/mrichter/AliRoot.git] / 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