]>
Commit | Line | Data |
---|---|---|
a3408ed3 | 1 | #include "TKDInterpolatorBase.h" |
2 | #include "TKDNodeInfo.h" | |
3 | #include "TKDTree.h" | |
4 | ||
b6f6b31a | 5 | #include "TROOT.h" |
a3408ed3 | 6 | #include "TClonesArray.h" |
7 | #include "TLinearFitter.h" | |
8 | #include "TTree.h" | |
9 | #include "TH2.h" | |
10 | #include "TObjArray.h" | |
11 | #include "TObjString.h" | |
a3408ed3 | 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 | //_________________________________________________________________ | |
1a439134 | 31 | TKDInterpolatorBase::TKDInterpolatorBase(Int_t dim) : |
50c4eb3a | 32 | fNSize(dim) |
afb669c1 | 33 | ,fNodes(NULL) |
34 | ,fNodesDraw(NULL) | |
b273c7cb | 35 | ,fStatus(0) |
50c4eb3a | 36 | ,fLambda(1 + dim + (dim*(dim+1)>>1)) |
37 | ,fDepth(-1) | |
38 | ,fAlpha(.5) | |
afb669c1 | 39 | ,fRefPoints(NULL) |
40 | ,fBuffer(NULL) | |
41 | ,fKDhelper(NULL) | |
42 | ,fFitter(NULL) | |
a3408ed3 | 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. | |
b273c7cb | 46 | UseWeights(); |
a3408ed3 | 47 | } |
48 | ||
49 | //_________________________________________________________________ | |
afb669c1 | 50 | Bool_t TKDInterpolatorBase::Build(Int_t n) |
a3408ed3 | 51 | { |
50c4eb3a | 52 | // allocate memory for data |
afb669c1 | 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 | } | |
a3408ed3 | 57 | |
afb669c1 | 58 | if(fNodes){ |
59 | Warning("TKDInterpolatorBase::Build()", "Data already allocated."); | |
60 | fNodes->Delete(); | |
61 | } else { | |
62 | fNodes = new TClonesArray("TKDNodeInfo", n); | |
63 | fNodes->SetOwner(); | |
53ee3cad | 64 | } |
afb669c1 | 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; | |
a3408ed3 | 90 | } |
91 | ||
92 | //_________________________________________________________________ | |
93 | TKDInterpolatorBase::~TKDInterpolatorBase() | |
94 | { | |
50c4eb3a | 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 | } | |
afb669c1 | 103 | if(fNodes){ |
104 | fNodes->Delete(); | |
105 | delete fNodes; | |
b6f6b31a | 106 | } |
afb669c1 | 107 | if(fNodesDraw) delete [] fNodesDraw; |
b6f6b31a | 108 | |
afb669c1 | 109 | TH2 *h2=NULL; |
110 | if((h2 = (TH2S*)gROOT->FindObject("hKDnodes"))) delete h2; | |
a3408ed3 | 111 | } |
112 | ||
113 | ||
114 | //__________________________________________________________________ | |
115 | Bool_t TKDInterpolatorBase::GetCOGPoint(Int_t inode, Float_t *&coord, Float_t &val, Float_t &err) const | |
116 | { | |
afb669c1 | 117 | if(inode < 0 || inode > GetNTNodes()) return kFALSE; |
a3408ed3 | 118 | |
afb669c1 | 119 | TKDNodeInfo *node = (TKDNodeInfo*)(*fNodes)[inode]; |
50c4eb3a | 120 | coord = &(node->Data()[0]); |
121 | val = node->Val()[0]; | |
122 | err = node->Val()[1]; | |
123 | return kTRUE; | |
a3408ed3 | 124 | } |
125 | ||
126 | //_________________________________________________________________ | |
1a439134 | 127 | TKDNodeInfo* TKDInterpolatorBase::GetNodeInfo(Int_t inode) const |
a3408ed3 | 128 | { |
afb669c1 | 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; | |
a3408ed3 | 137 | } |
138 | ||
b273c7cb | 139 | //_________________________________________________________________ |
140 | Bool_t TKDInterpolatorBase::GetRange(Int_t ax, Float_t &min, Float_t &max) const | |
141 | { | |
afb669c1 | 142 | if(!fNodes) return kFALSE; |
143 | Int_t ndim = ((TKDNodeInfo*)(*fNodes)[0])->GetDimension(); | |
b273c7cb | 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; | |
afb669c1 | 150 | for(Int_t in=GetNTNodes(); in--; ){ |
151 | TKDNodeInfo *node = (TKDNodeInfo*)((*fNodes)[in]); | |
b273c7cb | 152 | node->GetBoundary(ax, axmin, axmax); |
153 | if(axmin<min) min = axmin; | |
154 | if(axmax>max) max = axmax; | |
155 | } | |
156 | ||
157 | return kTRUE; | |
158 | } | |
a3408ed3 | 159 | |
160 | //__________________________________________________________________ | |
b273c7cb | 161 | void TKDInterpolatorBase::GetStatus(Option_t *opt) |
a3408ed3 | 162 | { |
163 | // Prints the status of the interpolator | |
164 | ||
b273c7cb | 165 | printf("Interpolator Status[%d] :\n", fStatus); |
50c4eb3a | 166 | printf(" Dim : %d [%d]\n", fNSize, fLambda); |
b273c7cb | 167 | printf(" Method : %s\n", UseCOG() ? "COG" : "INT"); |
168 | printf(" Store : %s\n", HasStore() ? "YES" : "NO"); | |
169 | printf(" Weights: %s\n", UseWeights() ? "YES" : "NO"); | |
50c4eb3a | 170 | |
b273c7cb | 171 | if(strcmp(opt, "all") != 0 ) return; |
afb669c1 | 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]; | |
50c4eb3a | 175 | printf("%d ", i); node->Print(); |
176 | } | |
a3408ed3 | 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) | |
afb669c1 | 188 | |
50c4eb3a | 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); | |
afb669c1 | 192 | if(nodeIndex<0){ |
193 | Error("TKDInterpolatorBase::Eval()", "Can not retrive node for data point."); | |
194 | result = 0.; | |
50c4eb3a | 195 | error = 1.E10; |
196 | return 0.; | |
197 | } | |
afb669c1 | 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 | } | |
50c4eb3a | 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++){ | |
afb669c1 | 209 | fRefPoints[id] = new Float_t[GetNTNodes()]; |
210 | for(int in=0; in<GetNTNodes(); in++) fRefPoints[id][in] = ((TKDNodeInfo*)(*fNodes)[in])->Data()[id]; | |
50c4eb3a | 211 | } |
afb669c1 | 212 | Info("TKDInterpolatorBase::Eval()", Form("Build TKDTree(%d, %d, %d)", GetNTNodes(), fNSize, kNhelper)); |
213 | fKDhelper = new TKDTreeIF(GetNTNodes(), fNSize, kNhelper, fRefPoints); | |
53ee3cad | 214 | fKDhelper->Build(); |
215 | fKDhelper->MakeBoundariesExact(); | |
50c4eb3a | 216 | } |
217 | if(!fFitter) fFitter = new TLinearFitter(fLambda, Form("hyp%d", fLambda-1)); | |
218 | ||
219 | // generate parabolic for nD | |
afb669c1 | 220 | //Float_t alpha = Float_t(2*lambda + 1) / GetNTNodes(); // the bandwidth or smoothing parameter |
221 | //Int_t npoints = Int_t(alpha * GetNTNodes()); | |
50c4eb3a | 222 | //printf("Params : %d NPoints %d\n", lambda, npoints); |
223 | // prepare workers | |
224 | ||
225 | Int_t ipar, // local looping variable | |
53ee3cad | 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 | |
50c4eb3a | 230 | d, // NN normalized distance |
231 | w0, // work | |
232 | w; // tri-cubic weight function | |
233 | ||
53ee3cad | 234 | Bool_t kDOWN = kFALSE; |
50c4eb3a | 235 | do{ |
b273c7cb | 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 | } | |
53ee3cad | 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; | |
20932c8c | 243 | // clean memory |
244 | delete [] dist; delete [] index; | |
53ee3cad | 245 | return 0.; |
246 | } else npoints = npoints_new; | |
afb669c1 | 247 | if(npoints > GetNTNodes()){ |
248 | Warning("TKDInterpolatorBase::Eval()", Form("The number of interpolation points requested (%d) exceeds number of PDF values (%d). Downscale.", npoints, GetNTNodes())); | |
249 | npoints = GetNTNodes(); | |
53ee3cad | 250 | kDOWN = kTRUE; |
251 | } | |
252 | ||
50c4eb3a | 253 | // find nearest neighbors |
254 | for(int idim=0; idim<fNSize; idim++) pointF[idim] = (Float_t)point[idim]; | |
255 | fKDhelper->FindNearestNeighbors(pointF, npoints+1, index, dist); | |
53ee3cad | 256 | |
50c4eb3a | 257 | // add points to fitter |
258 | fFitter->ClearPoints(); | |
afb669c1 | 259 | TKDNodeInfo *tnode = NULL; |
50c4eb3a | 260 | for(int in=0; in<npoints; in++){ |
afb669c1 | 261 | tnode = (TKDNodeInfo*)(*fNodes)[index[in]]; |
50c4eb3a | 262 | //tnode->Print(); |
b273c7cb | 263 | if(UseCOG()){ // COG |
264 | Float_t *p = &(tnode->Data()[0]); | |
265 | ipar = 0; | |
266 | for(int idim=0; idim<fNSize; idim++){ | |
267 | fBuffer[ipar++] = p[idim]; | |
268 | for(int jdim=idim; jdim<fNSize; jdim++) fBuffer[ipar++] = p[idim]*p[jdim]; | |
269 | } | |
270 | } else { // INT | |
50c4eb3a | 271 | Float_t *bounds = &(tnode->Data()[fNSize]); |
272 | ipar = 0; | |
273 | for(int idim=0; idim<fNSize; idim++){ | |
274 | fBuffer[ipar++] = .5*(bounds[2*idim] + bounds[2*idim+1]); | |
275 | fBuffer[ipar++] = (bounds[2*idim]*bounds[2*idim] + bounds[2*idim] * bounds[2*idim+1] + bounds[2*idim+1] * bounds[2*idim+1])/3.; | |
276 | 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; | |
277 | } | |
50c4eb3a | 278 | } |
279 | ||
280 | // calculate tri-cubic weighting function | |
b273c7cb | 281 | if(UseWeights()){ |
282 | d = dist[in]/dist[npoints]; | |
50c4eb3a | 283 | w0 = (1. - d*d*d); w = w0*w0*w0; |
b273c7cb | 284 | if(w<1.e-30) continue; |
50c4eb3a | 285 | } else w = 1.; |
286 | ||
b273c7cb | 287 | // printf("%2d d[%f] w[%f] x[", index[in], d, w); |
53ee3cad | 288 | // for(int idim=0; idim<fLambda-1; idim++) printf("%f ", fBuffer[idim]); |
b273c7cb | 289 | // printf("]\n"); printf("v[%f +- %f] (%f, %f)\n", tnode->Val()[0], tnode->Val()[1]/w, tnode->Val()[1], w); |
50c4eb3a | 290 | fFitter->AddPoint(fBuffer, tnode->Val()[0], tnode->Val()[1]/w); |
291 | } | |
53ee3cad | 292 | npoints_new = npoints+ (kDOWN ? 0 : kdN); |
50c4eb3a | 293 | } while(fFitter->Eval()); |
294 | delete [] index; | |
295 | delete [] dist; | |
296 | ||
297 | // retrive fitter results | |
298 | TMatrixD cov(fLambda, fLambda); | |
299 | TVectorD par(fLambda); | |
300 | fFitter->GetCovarianceMatrix(cov); | |
301 | fFitter->GetParameters(par); | |
302 | Double_t chi2 = fFitter->GetChisquare()/(npoints - 4 - fLambda); | |
303 | ||
304 | // store results | |
afb669c1 | 305 | node->Store(&par, HasStore()?&cov:NULL); |
50c4eb3a | 306 | |
307 | // Build df/dpi|x values | |
308 | Double_t *fdfdp = &fBuffer[fLambda]; | |
309 | ipar = 0; | |
310 | fdfdp[ipar++] = 1.; | |
311 | for(int idim=0; idim<fNSize; idim++){ | |
312 | fdfdp[ipar++] = point[idim]; | |
313 | for(int jdim=idim; jdim<fNSize; jdim++) fdfdp[ipar++] = point[idim]*point[jdim]; | |
314 | } | |
315 | ||
316 | // calculate estimation | |
317 | result =0.; error = 0.; | |
318 | for(int i=0; i<fLambda; i++){ | |
319 | result += fdfdp[i]*par(i); | |
320 | for(int j=0; j<fLambda; j++) error += fdfdp[i]*fdfdp[j]*cov(i,j); | |
321 | } | |
322 | error = TMath::Sqrt(error); | |
50c4eb3a | 323 | return chi2; |
a3408ed3 | 324 | } |
325 | ||
326 | //_________________________________________________________________ | |
b273c7cb | 327 | void TKDInterpolatorBase::DrawProjection(UInt_t ax1, UInt_t ax2) |
a3408ed3 | 328 | { |
329 | // Draw nodes structure projected on plane "ax1:ax2". The parameter | |
330 | // "depth" specifies the bucket size per node. If depth == -1 draw only | |
331 | // terminal nodes and evaluation points (default -1 i.e. bucket size per node equal bucket size specified by the user) | |
332 | // | |
50c4eb3a | 333 | |
b273c7cb | 334 | Float_t ax1min, ax1max, ax2min, ax2max; |
335 | GetRange(ax1, ax1min, ax1max); | |
336 | GetRange(ax2, ax2min, ax2max); | |
afb669c1 | 337 | TH2 *h2 = NULL; |
b273c7cb | 338 | if(!(h2 = (TH2S*)gROOT->FindObject("hKDnodes"))){ |
b6f6b31a | 339 | h2 = new TH2S("hKDnodes", "", 100, ax1min, ax1max, 100, ax2min, ax2max); |
b6f6b31a | 340 | } |
b273c7cb | 341 | h2->GetXaxis()->SetRangeUser(ax1min, ax1max); |
342 | h2->GetXaxis()->SetTitle(Form("x_{%d}", ax1)); | |
343 | h2->GetYaxis()->SetRangeUser(ax2min, ax2max); | |
344 | h2->GetYaxis()->SetTitle(Form("x_{%d}", ax2)); | |
50c4eb3a | 345 | h2->Draw(); |
b273c7cb | 346 | |
50c4eb3a | 347 | |
afb669c1 | 348 | if(!fNodesDraw) fNodesDraw = new TKDNodeInfo::TKDNodeDraw[GetNTNodes()]; |
349 | TKDNodeInfo::TKDNodeDraw *box = NULL; | |
350 | for(Int_t in=GetNTNodes(); in--; ){ | |
351 | box = &(fNodesDraw[in]); | |
352 | box->SetNode((TKDNodeInfo*)((*fNodes)[in]), fNSize, ax1, ax2); | |
b6f6b31a | 353 | box->Draw(); |
50c4eb3a | 354 | } |
b273c7cb | 355 | |
50c4eb3a | 356 | return; |
a3408ed3 | 357 | } |
358 | ||
53ee3cad | 359 | //_________________________________________________________________ |
360 | void TKDInterpolatorBase::SetAlpha(Float_t a) | |
361 | { | |
362 | if(a<0.5){ | |
363 | Warning("TKDInterpolatorBase::SetAlpha()", "The scale parameter has to be larger than 0.5"); | |
364 | fAlpha = 0.5; | |
365 | return; | |
366 | } | |
367 | // check value | |
afb669c1 | 368 | if(Int_t((a+1.)*fLambda) > GetNTNodes()){ |
369 | fAlpha = TMath::Max(0.5, Float_t(GetNTNodes())/fLambda-1.); | |
53ee3cad | 370 | Warning("TKDInterpolatorBase::SetAlpha()", Form("Interpolation neighborhood exceeds number of evaluation points. Downscale alpha to %f", fAlpha)); |
afb669c1 | 371 | printf("n[%d] nodes[%d]\n", Int_t((fAlpha+1.)*fLambda), GetNTNodes()); |
53ee3cad | 372 | return; |
373 | } | |
374 | fAlpha = a; | |
375 | return; | |
376 | } | |
377 |