]>
Commit | Line | Data |
---|---|---|
f2040a8f | 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 "TBox.h" | |
10 | #include "TGraph.h" | |
11 | #include "TMarker.h" | |
12 | ||
13 | ||
14 | ||
15 | ClassImp(TKDInterpolator) | |
16 | ||
17 | ///////////////////////////////////////////////////////////////////// | |
18 | // Memory setup of protected data memebers | |
19 | // fRefPoints : evaluation point of PDF for each terminal node of underlying KD Tree. | |
20 | // | 1st terminal node (fNDim point coordinates) | 2nd terminal node (fNDim point coordinates) | ... | |
21 | // | |
22 | // fRefValues : evaluation value/error of PDF for each terminal node of underlying KD Tree. | |
23 | // | 1st terminal node (value) | 2nd terminal node (value) | ... | 1st terminal node (error) | 2nd terminal node (error) | ... | |
24 | ///////////////////////////////////////////////////////////////////// | |
25 | ||
26 | //_________________________________________________________________ | |
27 | TKDInterpolator::TKDInterpolator() : TKDTreeIF() | |
28 | ,fNTNodes(0) | |
29 | ,fRefPoints(0x0) | |
30 | ,fRefValues(0x0) | |
31 | ,fDepth(-1) | |
32 | ,fTmpPoint(0x0) | |
33 | ,fKDhelper(0x0) | |
34 | ,fFitter(0x0) | |
35 | { | |
36 | } | |
37 | ||
38 | //_________________________________________________________________ | |
39 | TKDInterpolator::TKDInterpolator(Int_t npoints, Int_t ndim, UInt_t bsize, Float_t **data) : TKDTreeIF(npoints, ndim, bsize, data) | |
40 | ,fNTNodes(GetNTerminalNodes()) | |
41 | ,fRefPoints(0x0) | |
42 | ,fRefValues(0x0) | |
43 | ,fDepth(-1) | |
44 | ,fTmpPoint(0x0) | |
45 | ,fKDhelper(0x0) | |
46 | ,fFitter(0x0) | |
47 | { | |
48 | Build(); | |
49 | } | |
50 | ||
51 | ||
52 | //_________________________________________________________________ | |
53 | TKDInterpolator::TKDInterpolator(TTree *t, const Char_t *var, const Char_t *cut, UInt_t bsize) : TKDTreeIF() | |
54 | ,fNTNodes(0) | |
55 | ,fRefPoints(0x0) | |
56 | ,fRefValues(0x0) | |
57 | ,fDepth(-1) | |
58 | ,fTmpPoint(0x0) | |
59 | ,fKDhelper(0x0) | |
60 | ,fFitter(0x0) | |
61 | { | |
62 | // Alocate data from a tree. The variables which have to be analysed are | |
63 | // defined in the "var" parameter as a colon separated list. The format should | |
64 | // be identical to that used by TTree::Draw(). | |
65 | // | |
66 | // | |
67 | ||
68 | fNpoints = t->GetEntriesFast(); | |
69 | TObjArray *vars = TString(var).Tokenize(":"); | |
70 | fNDim = vars->GetEntriesFast(); | |
71 | if(fNDim > 6/*kDimMax*/) Warning("TKDInterpolator(TTree*, UInt_t, const Char_t)", Form("Variable number exceed maximum dimension %d. Results are unpredictable.", 6/*kDimMax*/)); | |
72 | fBucketSize = bsize; | |
73 | ||
74 | printf("Allocating %d points in %d dimensions.\n", fNpoints, fNDim); | |
75 | Float_t *mem = new Float_t[fNDim*fNpoints]; | |
76 | fData = new Float_t*[fNDim]; | |
77 | for(int idim=0; idim<fNDim; idim++) fData[idim] = &mem[idim*fNpoints]; | |
78 | kDataOwner = kTRUE; | |
79 | ||
80 | Double_t *v; | |
81 | for(int idim=0; idim<fNDim; idim++){ | |
82 | if(!(t->Draw(((TObjString*)(*vars)[idim])->GetName(), cut, "goff"))){ | |
83 | Warning("TKDInterpolator(TTree*, UInt_t, const Char_t)", Form("Can not access data for %s", ((TObjString*)(*vars)[idim])->GetName() )); | |
84 | continue; | |
85 | } | |
86 | v = t->GetV1(); | |
87 | for(int ip=0; ip<fNpoints; ip++) fData[idim][ip] = (Float_t)v[ip]; | |
88 | } | |
89 | TKDTreeIF::Build(); | |
90 | fNTNodes = GetNTerminalNodes(); | |
91 | Build(); | |
92 | } | |
93 | ||
94 | //_________________________________________________________________ | |
95 | TKDInterpolator::~TKDInterpolator() | |
96 | { | |
97 | if(fFitter) delete fFitter; | |
98 | if(fKDhelper) delete fKDhelper; | |
99 | if(fTmpPoint) delete [] fTmpPoint; | |
100 | ||
101 | if(fRefPoints){ | |
102 | for(int idim=0; idim<fNDim; idim++) delete [] fRefPoints[idim] ; | |
103 | delete [] fRefPoints; | |
104 | } | |
105 | if(fRefValues) delete [] fRefValues; | |
106 | } | |
107 | ||
108 | //_________________________________________________________________ | |
109 | void TKDInterpolator::Build() | |
110 | { | |
111 | if(!fBoundaries) MakeBoundaries(); | |
112 | ||
113 | // allocate memory for data | |
114 | fRefValues = new Float_t[fNTNodes]; | |
115 | fRefPoints = new Float_t*[fNDim]; | |
116 | for(int id=0; id<fNDim; id++){ | |
117 | fRefPoints[id] = new Float_t[fNTNodes]; | |
118 | for(int in=0; in<fNTNodes; in++) fRefPoints[id][in] = 0.; | |
119 | } | |
120 | ||
121 | Float_t *bounds = 0x0; | |
122 | Int_t *indexPoints; | |
123 | for(int inode=0, tnode = fNnodes; inode<fNTNodes-1; inode++, tnode++){ | |
124 | fRefValues[inode] = Float_t(fBucketSize)/fNpoints; | |
125 | bounds = GetBoundary(tnode); | |
126 | for(int idim=0; idim<fNDim; idim++) fRefValues[inode] /= (bounds[2*idim+1] - bounds[2*idim]); | |
127 | ||
128 | indexPoints = GetPointsIndexes(tnode); | |
129 | // loop points in this terminal node | |
130 | for(int idim=0; idim<fNDim; idim++){ | |
131 | for(int ip = 0; ip<fBucketSize; ip++) fRefPoints[idim][inode] += fData[idim][indexPoints[ip]]; | |
132 | fRefPoints[idim][inode] /= fBucketSize; | |
133 | } | |
134 | } | |
135 | ||
136 | // analyze last (incomplete) terminal node | |
137 | Int_t counts = fNpoints%fBucketSize; | |
138 | counts = counts ? counts : fBucketSize; | |
139 | Int_t inode = fNTNodes - 1, tnode = inode + fNnodes; | |
140 | fRefValues[inode] = Float_t(counts)/fNpoints; | |
141 | bounds = GetBoundary(tnode); | |
142 | for(int idim=0; idim<fNDim; idim++) fRefValues[inode] /= (bounds[2*idim+1] - bounds[2*idim]); | |
143 | ||
144 | indexPoints = GetPointsIndexes(tnode); | |
145 | // loop points in this terminal node | |
146 | for(int idim=0; idim<fNDim; idim++){ | |
147 | for(int ip = 0; ip<counts; ip++) fRefPoints[idim][inode] += fData[idim][indexPoints[ip]]; | |
148 | fRefPoints[idim][inode] /= counts; | |
149 | } | |
150 | } | |
151 | ||
152 | //_________________________________________________________________ | |
153 | Double_t TKDInterpolator::Eval(Float_t *point) | |
154 | { | |
155 | ||
156 | // calculate number of parameters in the parabolic expresion | |
157 | Int_t kNN = 1 + fNDim + fNDim*(fNDim+1)/2; | |
158 | ||
159 | // prepare workers | |
160 | if(!fTmpPoint) fTmpPoint = new Double_t[fNDim]; | |
161 | if(!fKDhelper) fKDhelper = new TKDTreeIF(GetNTerminalNodes(), fNDim, kNN, fRefPoints); | |
162 | if(!fFitter){ | |
163 | // generate formula for nD | |
164 | TString formula("1"); | |
165 | for(int idim=0; idim<fNDim; idim++){ | |
166 | formula += Form("++x[%d]", idim); | |
167 | for(int jdim=idim; jdim<fNDim; jdim++) formula += Form("++x[%d]*x[%d]", idim, jdim); | |
168 | } | |
169 | fFitter = new TLinearFitter(kNN, formula.Data()); | |
170 | } | |
171 | ||
172 | Int_t kNN_old = 0; | |
173 | Int_t *index; | |
174 | Float_t dist; | |
175 | fFitter->ClearPoints(); | |
176 | do{ | |
177 | if(!fKDhelper->FindNearestNeighbors(point, kNN, index, dist)){ | |
178 | Error("Eval()", Form("Failed retriving %d neighbours for point:", kNN)); | |
179 | for(int idim=0; idim<fNDim; idim++) printf("%f ", point[idim]); | |
180 | printf("\n"); | |
181 | return -1; | |
182 | } | |
183 | for(int in=kNN_old; in<kNN; in++){ | |
184 | for(int idim=0; idim<fNDim; idim++) fTmpPoint[idim] = fRefPoints[idim][index[in]]; | |
185 | fFitter->AddPoint(fTmpPoint, TMath::Log(fRefValues[index[in]]), 1.); | |
186 | } | |
187 | kNN_old = kNN; | |
188 | kNN += 4; | |
189 | } while(fFitter->Eval()); | |
190 | ||
191 | // calculate evaluation | |
192 | TVectorD par(kNN); | |
193 | fFitter->GetParameters(par); | |
194 | Double_t result = par[0]; | |
195 | Int_t ipar = 0; | |
196 | for(int idim=0; idim<fNDim; idim++){ | |
197 | result += par[++ipar]*point[idim]; | |
198 | for(int jdim=idim; jdim<fNDim; jdim++) result += par[++ipar]*point[idim]*point[jdim]; | |
199 | } | |
200 | return TMath::Exp(result); | |
201 | } | |
202 | ||
203 | ||
204 | //_________________________________________________________________ | |
205 | void TKDInterpolator::DrawNodes(Int_t depth, Int_t ax1, Int_t ax2) | |
206 | { | |
207 | // Draw nodes structure projected on plane "ax1:ax2". The parameter | |
208 | // "depth" specifies the bucket size per node. If depth == -1 draw only | |
209 | // terminal nodes and evaluation points (default -1 i.e. bucket size per node equal bucket size specified by the user) | |
210 | ||
211 | if(!fBoundaries) MakeBoundaries(); | |
212 | ||
213 | // Count nodes in specific view | |
214 | Int_t nnodes = 0; | |
215 | for(int inode = 0; inode <= 2*fNnodes; inode++){ | |
216 | if(depth == -1){ | |
217 | if(!IsTerminal(inode)) continue; | |
218 | } else if((inode+1) >> depth != 1) continue; | |
219 | nnodes++; | |
220 | } | |
221 | ||
222 | //printf("depth %d nodes %d\n", depth, nnodes); | |
223 | ||
224 | //TH2 *h2 = new TH2S("hframe", "", 100, fRange[2*ax1], fRange[2*ax1+1], 100, fRange[2*ax2], fRange[2*ax2+1]); | |
225 | TH2 *h2 = new TH2S("hframe", "", 100, 0., 1., 100, 0., 1.); | |
226 | h2->Draw(); | |
227 | ||
228 | const Float_t border = 0.;//1.E-4; | |
229 | TBox **node_array = new TBox*[nnodes], *node; | |
230 | Float_t *bounds = 0x0; | |
231 | nnodes = 0; | |
232 | for(int inode = 0; inode <= 2*fNnodes; inode++){ | |
233 | if(depth == -1){ | |
234 | if(!IsTerminal(inode)) continue; | |
235 | } else if((inode+1) >> depth != 1) continue; | |
236 | ||
237 | node = node_array[nnodes++]; | |
238 | bounds = GetBoundary(inode); | |
239 | node = new TBox(bounds[2*ax1]+border, bounds[2*ax2]+border, bounds[2*ax1+1]-border, bounds[2*ax2+1]-border); | |
240 | node->SetFillStyle(0); | |
241 | node->SetFillColor(51+inode); | |
242 | node->Draw(); | |
243 | } | |
244 | if(depth != -1) return; | |
245 | ||
246 | // Draw reference points | |
247 | TGraph *ref = new TGraph(GetNTerminalNodes()); | |
248 | ref->SetMarkerStyle(2); | |
249 | ref->SetMarkerColor(2); | |
250 | Float_t val, error; | |
251 | for(int inode = 0; inode < GetNTerminalNodes(); inode++) ref->SetPoint(inode, fRefPoints[ax1][inode], fRefPoints[ax2][inode]); | |
252 | ref->Draw("p"); | |
253 | return; | |
254 | } | |
255 | ||
256 | //_________________________________________________________________ | |
257 | void TKDInterpolator::DrawNode(Int_t tnode, Int_t ax1, Int_t ax2) | |
258 | { | |
259 | // Draw node "node" and the data points within. | |
260 | ||
261 | if(tnode < 0 || tnode >= GetNTerminalNodes()){ | |
262 | Warning("DrawNode()", Form("Terminal node %d outside defined range.", tnode)); | |
263 | return; | |
264 | } | |
265 | ||
266 | //TH2 *h2 = new TH2S("hframe", "", 1, fRange[2*ax1], fRange[2*ax1+1], 1, fRange[2*ax2], fRange[2*ax2+1]); | |
267 | TH2 *h2 = new TH2S("hframe", "", 1, 0., 1., 1, 0., 1.); | |
268 | h2->Draw(); | |
269 | ||
270 | Int_t inode = tnode; | |
271 | tnode += fNnodes; | |
272 | // select zone of interest in the indexes array | |
273 | Int_t *index = GetPointsIndexes(tnode); | |
274 | Int_t nPoints = (tnode == 2*fNnodes) ? fNpoints%fBucketSize : fBucketSize; | |
275 | ||
276 | printf("true index %d points %d\n", tnode, nPoints); | |
277 | ||
278 | // draw data points | |
279 | TGraph *g = new TGraph(nPoints); | |
280 | g->SetMarkerStyle(3); | |
281 | g->SetMarkerSize(.8); | |
282 | for(int ip = 0; ip<nPoints; ip++) g->SetPoint(ip, fData[ax1][index[ip]], fData[ax2][index[ip]]); | |
283 | g->Draw("p"); | |
284 | ||
285 | // draw estimation point | |
286 | TMarker *m=new TMarker(fRefPoints[ax1][inode], fRefPoints[ax2][inode], 2); | |
287 | m->SetMarkerColor(2); | |
288 | m->Draw(); | |
289 | ||
290 | // draw node contour | |
291 | Float_t *bounds = GetBoundary(tnode); | |
292 | TBox *n = new TBox(bounds[2*ax1], bounds[2*ax2], bounds[2*ax1+1], bounds[2*ax2+1]); | |
293 | n->SetFillStyle(0); | |
294 | n->Draw(); | |
295 | ||
296 | return; | |
297 | } | |
298 |