]>
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) |
33 | ,fNTNodes(0) | |
34 | ,fTNodes(0x0) | |
b6f6b31a | 35 | ,fTNodesDraw(0x0) |
b273c7cb | 36 | ,fStatus(0) |
50c4eb3a | 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) | |
a3408ed3 | 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. | |
b273c7cb | 47 | UseWeights(); |
a3408ed3 | 48 | } |
49 | ||
50 | //_________________________________________________________________ | |
1a439134 | 51 | void TKDInterpolatorBase::Build(Int_t n) |
a3408ed3 | 52 | { |
50c4eb3a | 53 | // allocate memory for data |
a3408ed3 | 54 | |
50c4eb3a | 55 | if(fTNodes) delete fTNodes; |
56 | fNTNodes = n; | |
53ee3cad | 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 | } | |
50c4eb3a | 61 | fTNodes = new TClonesArray("TKDNodeInfo", fNTNodes); |
62 | for(int in=0; in<fNTNodes; in++) new ((*fTNodes)[in]) TKDNodeInfo(fNSize); | |
a3408ed3 | 63 | } |
64 | ||
65 | //_________________________________________________________________ | |
66 | TKDInterpolatorBase::~TKDInterpolatorBase() | |
67 | { | |
50c4eb3a | 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 | } | |
b6f6b31a | 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; | |
a3408ed3 | 84 | } |
85 | ||
86 | ||
87 | //__________________________________________________________________ | |
88 | Bool_t TKDInterpolatorBase::GetCOGPoint(Int_t inode, Float_t *&coord, Float_t &val, Float_t &err) const | |
89 | { | |
50c4eb3a | 90 | if(inode < 0 || inode > fNTNodes) return kFALSE; |
a3408ed3 | 91 | |
50c4eb3a | 92 | TKDNodeInfo *node = (TKDNodeInfo*)(*fTNodes)[inode]; |
93 | coord = &(node->Data()[0]); | |
94 | val = node->Val()[0]; | |
95 | err = node->Val()[1]; | |
96 | return kTRUE; | |
a3408ed3 | 97 | } |
98 | ||
99 | //_________________________________________________________________ | |
1a439134 | 100 | TKDNodeInfo* TKDInterpolatorBase::GetNodeInfo(Int_t inode) const |
a3408ed3 | 101 | { |
50c4eb3a | 102 | if(!fTNodes || inode >= fNTNodes) return 0x0; |
103 | return (TKDNodeInfo*)(*fTNodes)[inode]; | |
a3408ed3 | 104 | } |
105 | ||
b273c7cb | 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 | } | |
a3408ed3 | 126 | |
127 | //__________________________________________________________________ | |
b273c7cb | 128 | void TKDInterpolatorBase::GetStatus(Option_t *opt) |
a3408ed3 | 129 | { |
130 | // Prints the status of the interpolator | |
131 | ||
b273c7cb | 132 | printf("Interpolator Status[%d] :\n", fStatus); |
50c4eb3a | 133 | printf(" Dim : %d [%d]\n", fNSize, fLambda); |
b273c7cb | 134 | printf(" Method : %s\n", UseCOG() ? "COG" : "INT"); |
135 | printf(" Store : %s\n", HasStore() ? "YES" : "NO"); | |
136 | printf(" Weights: %s\n", UseWeights() ? "YES" : "NO"); | |
50c4eb3a | 137 | |
b273c7cb | 138 | if(strcmp(opt, "all") != 0 ) return; |
50c4eb3a | 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 | } | |
a3408ed3 | 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) | |
50c4eb3a | 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){ | |
53ee3cad | 160 | Error("TKDInterpolatorBase::Eval()", "Can not retrive node for data point."); |
50c4eb3a | 161 | result = 0.; |
162 | error = 1.E10; | |
163 | return 0.; | |
164 | } | |
165 | TKDNodeInfo *node = (TKDNodeInfo*)(*fTNodes)[nodeIndex]; | |
b273c7cb | 166 | if(node->Cov() && !force) return node->CookPDF(point, result, error); |
50c4eb3a | 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 | } | |
b273c7cb | 176 | Info("TKDInterpolatorBase::Eval()", Form("Build TKDTree(%d, %d, %d)", fNTNodes, fNSize, kNhelper)); |
53ee3cad | 177 | fKDhelper = new TKDTreeIF(fNTNodes, fNSize, kNhelper, fRefPoints); |
178 | fKDhelper->Build(); | |
179 | fKDhelper->MakeBoundariesExact(); | |
50c4eb3a | 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 | |
53ee3cad | 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 | |
50c4eb3a | 194 | d, // NN normalized distance |
195 | w0, // work | |
196 | w; // tri-cubic weight function | |
197 | ||
53ee3cad | 198 | Bool_t kDOWN = kFALSE; |
50c4eb3a | 199 | do{ |
b273c7cb | 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 | } | |
53ee3cad | 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 | ||
50c4eb3a | 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); | |
53ee3cad | 218 | |
50c4eb3a | 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(); | |
b273c7cb | 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 | |
50c4eb3a | 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 | } | |
50c4eb3a | 240 | } |
241 | ||
242 | // calculate tri-cubic weighting function | |
b273c7cb | 243 | if(UseWeights()){ |
244 | d = dist[in]/dist[npoints]; | |
50c4eb3a | 245 | w0 = (1. - d*d*d); w = w0*w0*w0; |
b273c7cb | 246 | if(w<1.e-30) continue; |
50c4eb3a | 247 | } else w = 1.; |
248 | ||
b273c7cb | 249 | // printf("%2d d[%f] w[%f] x[", index[in], d, w); |
53ee3cad | 250 | // for(int idim=0; idim<fLambda-1; idim++) printf("%f ", fBuffer[idim]); |
b273c7cb | 251 | // printf("]\n"); printf("v[%f +- %f] (%f, %f)\n", tnode->Val()[0], tnode->Val()[1]/w, tnode->Val()[1], w); |
50c4eb3a | 252 | fFitter->AddPoint(fBuffer, tnode->Val()[0], tnode->Val()[1]/w); |
253 | } | |
53ee3cad | 254 | npoints_new = npoints+ (kDOWN ? 0 : kdN); |
50c4eb3a | 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 | |
b273c7cb | 267 | if(HasStore()) node->Store(par, cov); |
50c4eb3a | 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; | |
a3408ed3 | 287 | } |
288 | ||
289 | //_________________________________________________________________ | |
b273c7cb | 290 | void TKDInterpolatorBase::DrawProjection(UInt_t ax1, UInt_t ax2) |
a3408ed3 | 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 | // | |
50c4eb3a | 296 | |
b273c7cb | 297 | Float_t ax1min, ax1max, ax2min, ax2max; |
298 | GetRange(ax1, ax1min, ax1max); | |
299 | GetRange(ax2, ax2min, ax2max); | |
b6f6b31a | 300 | TH2 *h2 = 0x0; |
b273c7cb | 301 | if(!(h2 = (TH2S*)gROOT->FindObject("hKDnodes"))){ |
b6f6b31a | 302 | h2 = new TH2S("hKDnodes", "", 100, ax1min, ax1max, 100, ax2min, ax2max); |
b6f6b31a | 303 | } |
b273c7cb | 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)); | |
50c4eb3a | 308 | h2->Draw(); |
b273c7cb | 309 | |
50c4eb3a | 310 | |
b6f6b31a | 311 | if(!fTNodesDraw) fTNodesDraw = new TKDNodeInfo::TKDNodeDraw[fNTNodes]; |
312 | TKDNodeInfo::TKDNodeDraw *box = 0x0; | |
b273c7cb | 313 | for(Int_t in=fNTNodes; in--; ){ |
b6f6b31a | 314 | box = &(fTNodesDraw[in]); |
b273c7cb | 315 | box->SetNode((TKDNodeInfo*)((*fTNodes)[in]), fNSize, ax1, ax2); |
b6f6b31a | 316 | box->Draw(); |
50c4eb3a | 317 | } |
b273c7cb | 318 | |
50c4eb3a | 319 | return; |
a3408ed3 | 320 | } |
321 | ||
53ee3cad | 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 |