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