1 #include "TKDInterpolator.h"
3 #include "TLinearFitter.h"
8 #include "TObjString.h"
18 ClassImp(TKDInterpolator)
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) | ...
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 /////////////////////////////////////////////////////////////////////
29 //_________________________________________________________________
30 TKDInterpolator::TKDInterpolator() : TKDTreeIF()
39 // Default constructor. To be used with care since in this case building
40 // of data structure is completly left to the user responsability.
43 //_________________________________________________________________
44 TKDInterpolator::TKDInterpolator(Int_t npoints, Int_t ndim, UInt_t bsize, Float_t **data) : TKDTreeIF(npoints, ndim, bsize, data)
45 ,fNTNodes(GetNTerminalNodes())
53 // Wrapper constructor for the similar TKDTree one.
59 //_________________________________________________________________
60 TKDInterpolator::TKDInterpolator(TTree *t, const Char_t *var, const Char_t *cut, UInt_t bsize) : TKDTreeIF()
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().
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*/));
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() ));
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];
96 for(int ip=0; ip<fNpoints; ip++) fData[idim][ip] = (Float_t)v[ip];
99 fNTNodes = GetNTerminalNodes();
103 //_________________________________________________________________
104 TKDInterpolator::~TKDInterpolator()
106 if(fFitter) delete fFitter;
107 if(fKDhelper) delete fKDhelper;
108 if(fTmpPoint) delete [] fTmpPoint;
111 for(int idim=0; idim<fNDim; idim++) delete [] fRefPoints[idim] ;
112 delete [] fRefPoints;
114 if(fRefValues) delete [] fRefValues;
117 //_________________________________________________________________
118 void TKDInterpolator::Build()
120 // Fill interpolator's data array i.e.
121 // - estimation points
122 // - corresponding PDF values
124 if(!fBoundaries) MakeBoundaries();
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.;
134 Float_t *bounds = 0x0;
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]);
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;
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]);
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;
165 //_________________________________________________________________
166 Double_t TKDInterpolator::Eval(const Double_t *point, Int_t npoints)
168 // Evaluate PDF at k-dimensional position "point". The initial number of
169 // neighbour estimation points is set to "npoints"
171 //Int_t npoints = Int_t(alpha * fNTNodes);
172 //printf("Params : %d NPoints %d\n", lambda, npoints);
174 if(!fTmpPoint) fTmpPoint = new Double_t[fNDim];
175 if(!fKDhelper) fKDhelper = new TKDTreeIF(GetNTerminalNodes(), fNDim, npoints, fRefPoints);
177 // generate parabolic for nD
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);
187 fFitter = new TLinearFitter(lambda, formula.Data());
188 Info("Eval(const Double_t*, Int_t)", Form("Using %s for local interpolation.", formula.Data()));
192 for(int idim=0; idim<fNDim; idim++) pointF[idim] = point[idim];
195 Float_t dist, d0, w0, w;
196 Double_t uncertainty = TMath::Sqrt(1./fBucketSize);
197 fFitter->ClearPoints();
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]);
205 for(int in=istart; in<npoints; in++){
206 //printf("%d index[%2d] x(", in, index[in]);
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]);
214 w0 = (1. - d0*d0*d0);
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);
222 } while(fFitter->Eval());
224 // calculate evaluation
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];
231 //printf("\tResult : %f [%f]\n", TMath::Exp(result), result);
232 return TMath::Exp(result);
236 //_________________________________________________________________
237 void TKDInterpolator::DrawNodes(UInt_t ax1, UInt_t ax2, Int_t depth)
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)
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 !
248 if(!fBoundaries) MakeBoundaries();
250 // Count nodes in specific view
252 for(int inode = 0; inode <= 2*fNnodes; inode++){
254 if(!IsTerminal(inode)) continue;
255 } else if((inode+1) >> depth != 1) continue;
259 //printf("depth %d nodes %d\n", depth, nnodes);
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));
267 const Float_t border = 0.;//1.E-4;
268 TBox *node_array = new TBox[nnodes], *node;
269 Float_t *bounds = 0x0;
271 for(int inode = 0; inode <= 2*fNnodes; inode++){
273 if(!IsTerminal(inode)) continue;
274 } else if((inode+1) >> depth != 1) continue;
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);
283 if(depth != -1) return;
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]);
295 //_________________________________________________________________
296 void TKDInterpolator::DrawNode(Int_t tnode, UInt_t ax1, UInt_t ax2)
298 // Draw node "node" and the data points within.
301 // This function creates some graphical objects
302 // but don't delete it. Abusing this function may cause memory leaks !
304 if(tnode < 0 || tnode >= GetNTerminalNodes()){
305 Warning("DrawNode()", Form("Terminal node %d outside defined range.", tnode));
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;
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]]);
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);
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]);
330 if(gPad) gPad->Clear();