alocate memory dynamically not on stack
[u/mrichter/AliRoot.git] / STAT / TKDInterpolator.cxx
1 #include "TKDInterpolator.h"
2
3 #include "TLinearFitter.h"
4 #include "TVector.h"
5 #include "TTree.h"
6 #include "TH2.h"
7 #include "TObjArray.h"
8 #include "TObjString.h"
9 #include "TPad.h"
10 #include "TBox.h"
11 #include "TGraph.h"
12 #include "TMarker.h"
13 #include "TRandom.h"
14 #include "TROOT.h"
15
16
17
18 ClassImp(TKDInterpolator)
19
20 /////////////////////////////////////////////////////////////////////
21 // Memory setup of protected data memebers
22 // fRefPoints : evaluation point of PDF for each terminal node of underlying KD Tree.
23 // | 1st terminal node (fNDim point coordinates) | 2nd terminal node (fNDim point coordinates) | ...
24 //
25 // fRefValues : evaluation value/error of PDF for each terminal node of underlying KD Tree.
26 // | 1st terminal node (value) | 2nd terminal node (value) | ... | 1st terminal node (error) | 2nd terminal node (error) | ...
27 /////////////////////////////////////////////////////////////////////
28
29 //_________________________________________________________________
30 TKDInterpolator::TKDInterpolator() : TKDTreeIF()
31         ,fNTNodes(0)
32         ,fRefPoints(0x0)
33         ,fRefValues(0x0)
34         ,fDepth(-1)
35         ,fTmpPoint(0x0)
36         ,fKDhelper(0x0)
37         ,fFitter(0x0)
38 {
39 // Default constructor. To be used with care since in this case building
40 // of data structure is completly left to the user responsability.
41 }
42
43 //_________________________________________________________________
44 TKDInterpolator::TKDInterpolator(Int_t npoints, Int_t ndim, UInt_t bsize, Float_t **data) : TKDTreeIF(npoints, ndim, bsize, data)
45         ,fNTNodes(GetNTerminalNodes())
46         ,fRefPoints(0x0)
47         ,fRefValues(0x0)
48         ,fDepth(-1)
49         ,fTmpPoint(0x0)
50         ,fKDhelper(0x0)
51         ,fFitter(0x0)
52 {
53 // Wrapper constructor for the similar TKDTree one.
54         
55         Build();
56 }
57
58
59 //_________________________________________________________________
60 TKDInterpolator::TKDInterpolator(TTree *t, const Char_t *var, const Char_t *cut, UInt_t bsize) : TKDTreeIF()
61         ,fNTNodes(0)
62         ,fRefPoints(0x0)
63         ,fRefValues(0x0)
64         ,fDepth(-1)
65         ,fTmpPoint(0x0)
66         ,fKDhelper(0x0)
67         ,fFitter(0x0)
68 {
69 // Alocate data from a tree. The variables which have to be analysed are
70 // defined in the "var" parameter as a colon separated list. The format should
71 // be identical to that used by TTree::Draw().
72 //
73 // 
74
75         TObjArray *vars = TString(var).Tokenize(":");
76         fNDim = vars->GetEntriesFast();
77         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*/));
78         fBucketSize = bsize;
79
80         Int_t np;
81         Double_t *v;
82         for(int idim=0; idim<fNDim; idim++){
83                 if(!(np = t->Draw(((TObjString*)(*vars)[idim])->GetName(), cut, "goff"))){
84                         Warning("TKDInterpolator(TTree*, const Char_t, const Char_t, UInt_t)", Form("Can not access data for %s", ((TObjString*)(*vars)[idim])->GetName() ));
85                         continue;
86                 }
87                 if(!fNpoints){
88                         fNpoints = np;
89                         Info("TKDInterpolator(TTree*, const Char_t, const Char_t, UInt_t)", Form("Allocating %d data points in %d dimensions.", fNpoints, fNDim));
90                         //Float_t *mem = new Float_t[fNDim*fNpoints];
91                         fData = new Float_t*[fNDim];
92                         for(int idim=0; idim<fNDim; idim++) fData[idim] = new Float_t[fNpoints]; //&mem[idim*fNpoints];
93                         kDataOwner = kTRUE;
94                 }
95                 v = t->GetV1();
96                 for(int ip=0; ip<fNpoints; ip++) fData[idim][ip] = (Float_t)v[ip];
97         }
98         TKDTreeIF::Build();
99         fNTNodes = GetNTerminalNodes();
100         Build();
101 }
102
103 //_________________________________________________________________
104 TKDInterpolator::~TKDInterpolator()
105 {
106         if(fFitter) delete fFitter;
107         if(fKDhelper) delete fKDhelper;
108         if(fTmpPoint) delete [] fTmpPoint;
109         
110         if(fRefPoints){
111                 for(int idim=0; idim<fNDim; idim++) delete [] fRefPoints[idim] ;
112                 delete [] fRefPoints;
113         }
114         if(fRefValues) delete [] fRefValues;
115 }
116
117 //_________________________________________________________________
118 void TKDInterpolator::Build()
119 {
120 // Fill interpolator's data array i.e.
121 //  - estimation points 
122 //  - corresponding PDF values
123
124         if(!fBoundaries) MakeBoundaries();
125         
126         // allocate memory for data
127         fRefValues = new Float_t[fNTNodes];
128         fRefPoints = new Float_t*[fNDim];
129         for(int id=0; id<fNDim; id++){
130                 fRefPoints[id] = new Float_t[fNTNodes];
131                 for(int in=0; in<fNTNodes; in++) fRefPoints[id][in] = 0.;
132         }
133
134         Float_t *bounds = 0x0;
135         Int_t *indexPoints;
136         for(int inode=0, tnode = fNnodes; inode<fNTNodes-1; inode++, tnode++){
137                 fRefValues[inode] =  Float_t(fBucketSize)/fNpoints;
138                 bounds = GetBoundary(tnode);
139                 for(int idim=0; idim<fNDim; idim++) fRefValues[inode] /= (bounds[2*idim+1] - bounds[2*idim]);
140
141                 indexPoints = GetPointsIndexes(tnode);
142                 // loop points in this terminal node
143                 for(int idim=0; idim<fNDim; idim++){
144                         for(int ip = 0; ip<fBucketSize; ip++) fRefPoints[idim][inode] += fData[idim][indexPoints[ip]];
145                         fRefPoints[idim][inode] /= fBucketSize;
146                 }
147         }
148
149         // analyze last (incomplete) terminal node
150         Int_t counts = fNpoints%fBucketSize;
151         counts = counts ? counts : fBucketSize;
152         Int_t inode = fNTNodes - 1, tnode = inode + fNnodes;
153         fRefValues[inode] =  Float_t(counts)/fNpoints;
154         bounds = GetBoundary(tnode);
155         for(int idim=0; idim<fNDim; idim++) fRefValues[inode] /= (bounds[2*idim+1] - bounds[2*idim]);
156
157         indexPoints = GetPointsIndexes(tnode);
158         // loop points in this terminal node
159         for(int idim=0; idim<fNDim; idim++){
160                 for(int ip = 0; ip<counts; ip++) fRefPoints[idim][inode] += fData[idim][indexPoints[ip]];
161                 fRefPoints[idim][inode] /= counts;
162         }
163 }
164
165 //_________________________________________________________________
166 Double_t TKDInterpolator::Eval(const Double_t *point, Int_t npoints)
167 {
168 // Evaluate PDF at k-dimensional position "point". The initial number of
169 // neighbour estimation points is set to "npoints"
170         
171         //Int_t npoints = Int_t(alpha * fNTNodes);
172         //printf("Params : %d NPoints %d\n", lambda, npoints);
173         // prepare workers
174         if(!fTmpPoint) fTmpPoint = new Double_t[fNDim];
175         if(!fKDhelper) fKDhelper = new TKDTreeIF(GetNTerminalNodes(), fNDim, npoints, fRefPoints);
176         if(!fFitter){
177                 // generate parabolic for nD
178                 
179                 // calculate number of parameters in the parabolic expresion
180                 Int_t lambda = 1 + fNDim + fNDim*(fNDim+1)/2;
181                 //Float_t alpha = Float_t(2*lambda + 1) / fNTNodes; // the bandwidth or smoothing parameter
182                 TString formula("1");
183                 for(int idim=0; idim<fNDim; idim++){
184                         formula += Form("++x[%d]", idim);
185                         for(int jdim=idim; jdim<fNDim; jdim++) formula += Form("++x[%d]*x[%d]", idim, jdim);
186                 }
187                 fFitter = new TLinearFitter(lambda, formula.Data());
188                 Info("Eval(const Double_t*, Int_t)", Form("Using %s for local interpolation.", formula.Data()));
189         }
190
191         Float_t pointF[50];
192         for(int idim=0; idim<fNDim; idim++) pointF[idim] = point[idim];
193         Int_t istart = 0;
194         Int_t *index;
195         Float_t dist, d0, w0, w;
196         Double_t uncertainty = TMath::Sqrt(1./fBucketSize); 
197         fFitter->ClearPoints();
198         do{
199                 if(!fKDhelper->FindNearestNeighbors(pointF, npoints+1, index, dist)){
200                         Error("Eval()", Form("Failed retriving %d neighbours for point:", npoints));
201                         for(int idim=0; idim<fNDim; idim++) printf("%f ", point[idim]);
202                         printf("\n");
203                         return -1;
204                 }
205                 for(int in=istart; in<npoints; in++){
206                         //printf("%d index[%2d] x(", in, index[in]);
207                         d0 = 0.;
208                         for(int idim=0; idim<fNDim; idim++){
209                                 fTmpPoint[idim] = fRefPoints[idim][index[in]];
210                                 //printf("%6.4f ", fTmpPoint[idim]);
211                                 d0 += TMath::Abs(fTmpPoint[idim] - point[idim]);
212                         }
213                         d0 /= dist;
214                         w0 = (1. - d0*d0*d0);
215                         w = w0*w0*w0;
216
217                         //printf(") f = %f [%f] d0 = %6.4f   w = %6.4f  \n", fRefValues[index[in]], TMath::Log(fRefValues[index[in]]), d0, w);
218                         fFitter->AddPoint(fTmpPoint, TMath::Log(fRefValues[index[in]]), uncertainty/w);
219                 }
220                 istart = npoints;
221                 npoints += 4;
222         } while(fFitter->Eval());
223
224         // calculate evaluation
225         Int_t ipar = 0;
226         Double_t result = fFitter->GetParameter(ipar++);
227         for(int idim=0; idim<fNDim; idim++){
228                 result += fFitter->GetParameter(ipar++)*point[idim];
229                 for(int jdim=idim; jdim<fNDim; jdim++) result += fFitter->GetParameter(ipar++)*point[idim]*point[jdim];
230         }
231         //printf("\tResult : %f [%f]\n", TMath::Exp(result), result);
232         return TMath::Exp(result);
233 }
234
235
236 //_________________________________________________________________
237 void TKDInterpolator::DrawNodes(UInt_t ax1, UInt_t ax2, Int_t depth)
238 {
239 // Draw nodes structure projected on plane "ax1:ax2". The parameter
240 // "depth" specifies the bucket size per node. If depth == -1 draw only
241 // terminal nodes and evaluation points (default -1 i.e. bucket size per node equal bucket size specified by the user)
242 //
243 // Observation:
244 // This function creates the nodes (TBox) array for the specified depth
245 // but don't delete it. Abusing this function may cause memory leaks !
246
247
248         if(!fBoundaries) MakeBoundaries();
249
250         // Count nodes in specific view
251         Int_t nnodes = 0;
252         for(int inode = 0; inode <= 2*fNnodes; inode++){
253                 if(depth == -1){
254                         if(!IsTerminal(inode)) continue;
255                 } else if((inode+1) >> depth != 1) continue;
256                 nnodes++;
257         }
258
259         //printf("depth %d nodes %d\n", depth, nnodes);
260         
261         TH2 *h2 = 0x0;
262         if(!(h2 = (TH2S*)gROOT->FindObject("hNodes"))) h2 = new TH2S("hNodes", "", 100, fRange[2*ax1], fRange[2*ax1+1], 100, fRange[2*ax2], fRange[2*ax2+1]);
263         h2->GetXaxis()->SetTitle(Form("x_{%d}", ax1));
264         h2->GetYaxis()->SetTitle(Form("x_{%d}", ax2));
265         h2->Draw();
266         
267         const Float_t border = 0.;//1.E-4;
268         TBox *node_array = new TBox[nnodes], *node;
269         Float_t *bounds = 0x0;
270         nnodes = 0;
271         for(int inode = 0; inode <= 2*fNnodes; inode++){
272                 if(depth == -1){
273                         if(!IsTerminal(inode)) continue;
274                 } else if((inode+1) >> depth != 1) continue;
275
276                 node = &node_array[nnodes++];
277                 //node = new TBox(bounds[2*ax1]+border, bounds[2*ax2]+border, bounds[2*ax1+1]-border, bounds[2*ax2+1]-border);
278                 node->SetFillStyle(3002);       
279                 node->SetFillColor(50+Int_t(gRandom->Uniform()*50.));
280                 bounds = GetBoundary(inode);
281                 node->DrawBox(bounds[2*ax1]+border, bounds[2*ax2]+border, bounds[2*ax1+1]-border, bounds[2*ax2+1]-border);
282         }
283         if(depth != -1) return;
284
285         // Draw reference points
286         TGraph *ref = new TGraph(GetNTerminalNodes());
287         ref->SetMarkerStyle(3);
288         ref->SetMarkerSize(.7);
289         ref->SetMarkerColor(2);
290         for(int inode = 0; inode < GetNTerminalNodes(); inode++) ref->SetPoint(inode, fRefPoints[ax1][inode], fRefPoints[ax2][inode]);
291         ref->Draw("p");
292         return;
293 }
294
295 //_________________________________________________________________
296 void TKDInterpolator::DrawNode(Int_t tnode, UInt_t ax1, UInt_t ax2)
297 {
298 // Draw node "node" and the data points within.
299 //
300 // Observation:
301 // This function creates some graphical objects
302 // but don't delete it. Abusing this function may cause memory leaks !
303
304         if(tnode < 0 || tnode >= GetNTerminalNodes()){
305                 Warning("DrawNode()", Form("Terminal node %d outside defined range.", tnode));
306                 return;
307         }
308
309         Int_t inode = tnode;
310         tnode += fNnodes;
311         // select zone of interest in the indexes array
312         Int_t *index = GetPointsIndexes(tnode);
313         Int_t nPoints = (tnode == 2*fNnodes) ? fNpoints%fBucketSize : fBucketSize;
314
315         // draw data points
316         TGraph *g = new TGraph(nPoints);
317         g->SetMarkerStyle(7);
318         for(int ip = 0; ip<nPoints; ip++) g->SetPoint(ip, fData[ax1][index[ip]], fData[ax2][index[ip]]);
319
320         // draw estimation point
321         TMarker *m=new TMarker(fRefPoints[ax1][inode], fRefPoints[ax2][inode], 20);
322         m->SetMarkerColor(2);
323         m->SetMarkerSize(1.7);
324         
325         // draw node contour
326         Float_t *bounds = GetBoundary(tnode);
327         TBox *n = new TBox(bounds[2*ax1], bounds[2*ax2], bounds[2*ax1+1], bounds[2*ax2+1]);
328         n->SetFillStyle(0);
329
330         if(gPad) gPad->Clear(); 
331         g->Draw("ap");
332         m->Draw();
333         n->Draw();
334         
335         return;
336 }
337