]> git.uio.no Git - u/mrichter/AliRoot.git/blame - STAT/TKDTree.cxx
git://git.uio.no / u / mrichter / AliRoot.git / blame
? search:
summary | shortlog | log | commit | commitdiff | tree
blob | blame (incremental) | history | HEAD
New protection against noise clusters in not existing modules + bug fix in geometry...
[u/mrichter/AliRoot.git] / STAT / TKDTree.cxx
CommitLineData
f2040a8f 1#include "TKDTree.h"
2
3#include "TString.h"
4#include <string.h>
5
6#ifndef R__ALPHA
7templateClassImp(TKDTree)
8#endif
9//////////////////////////////////////////////////////////////////////
10// Memory setup of protected data members:
11//
12// kDataOwner; // Toggle ownership of the data
13// fNnodes:
14// size of branch node array, and index of first terminal node
15// fNDim; // number of variables
16// fNpoints; // number of multidimensional points
17// fBucketSize; // number of data points per terminal node
18//
19// fIndPoints : array containing rearanged indexes according to nodes:
20// | point indexes from 1st terminal node (fBucketSize elements) | ... | point indexes from the last terminal node (<= fBucketSize elements)
21//
22// Value **fData;
23// fRange : array containing the boundaries of the domain:
24// | 1st dimension (min + max) | 2nd dimension (min + max) | ...
25// fBoundaries : nodes boundaries
26// | 1st node {1st dim * 2 elements | 2nd dim * 2 elements | ...} | 2nd node {...} | ...
27// the nodes are arranged in the following order :
28// - first fNnodes are primary nodes
29// - after are the terminal nodes
30// fNodes : array of primary nodes
31///////////////////////////////////////////////////////////////////////
117c62ab 32
117c62ab 33
f2040a8f 34//_________________________________________________________________
35template <typename Index, typename Value>
36TKDTree<Index, Value>::TKDTree() :
37 TObject()
38 ,kDataOwner(kFALSE)
39 ,fNnodes(0)
40 ,fNDim(0)
316a7f5a 41 ,fNDimm(0)
f2040a8f 42 ,fNpoints(0)
43 ,fBucketSize(0)
117c62ab 44 ,fAxis(0x0)
45 ,fValue(0x0)
f2040a8f 46 ,fRange(0x0)
316a7f5a 47 ,fData(0x0)
f2040a8f 48 ,fBoundaries(0x0)
316a7f5a 49 ,fkNNdim(0)
f2040a8f 50 ,fkNN(0x0)
316a7f5a 51 ,fkNNdist(0x0)
52 ,fDistBuffer(0x0)
53 ,fIndBuffer(0x0)
a9c20b1f 54 ,fIndPoints(0x0)
55 ,fRowT0(0)
56 ,fCrossNode(0)
57 ,fOffset(0)
f2040a8f 58{
59// Default constructor. Nothing is built
60}
61
62
63//_________________________________________________________________
64template <typename Index, typename Value>
65TKDTree<Index, Value>::TKDTree(Index npoints, Index ndim, UInt_t bsize, Value **data) :
66 TObject()
67 ,kDataOwner(kFALSE)
68 ,fNnodes(0)
69 ,fNDim(ndim)
316a7f5a 70 ,fNDimm(2*ndim)
f2040a8f 71 ,fNpoints(npoints)
72 ,fBucketSize(bsize)
117c62ab 73 ,fAxis(0x0)
74 ,fValue(0x0)
f2040a8f 75 ,fRange(0x0)
316a7f5a 76 ,fData(data) //Columnwise!!!!!
f2040a8f 77 ,fBoundaries(0x0)
316a7f5a 78 ,fkNNdim(0)
f2040a8f 79 ,fkNN(0x0)
316a7f5a 80 ,fkNNdist(0x0)
81 ,fDistBuffer(0x0)
82 ,fIndBuffer(0x0)
a9c20b1f 83 ,fIndPoints(0x0)
84 ,fRowT0(0)
85 ,fCrossNode(0)
86 ,fOffset(0)
f2040a8f 87{
316a7f5a 88// Allocate data by hand. See TKDTree(TTree*, const Char_t*) for an automatic method.
89
f2040a8f 90 Build();
91}
92
f2040a8f 93//_________________________________________________________________
94template <typename Index, typename Value>
95TKDTree<Index, Value>::~TKDTree()
96{
316a7f5a 97 if (fIndBuffer) delete [] fIndBuffer;
98 if (fDistBuffer) delete [] fDistBuffer;
99 if (fkNNdist) delete [] fkNNdist;
f2040a8f 100 if (fkNN) delete [] fkNN;
117c62ab 101 if (fAxis) delete [] fAxis;
102 if (fValue) delete [] fValue;
f2040a8f 103 if (fIndPoints) delete [] fIndPoints;
104 if (fRange) delete [] fRange;
105 if (fBoundaries) delete [] fBoundaries;
106 if (kDataOwner && fData){
107 for(int idim=0; idim<fNDim; idim++) delete [] fData[idim];
108 delete [] fData;
109 }
110}
111
112
113//_________________________________________________________________
114template <typename Index, typename Value>
115void TKDTree<Index, Value>::Build(){
116 //
117 // Build binning
118 //
119 // 1. calculate number of nodes
120 // 2. calculate first terminal row
121 // 3. initialize index array
122 // 4. non recursive building of the binary tree
123 //
124 //1.
125 fNnodes = fNpoints/fBucketSize-1;
126 if (fNpoints%fBucketSize) fNnodes++;
127
128 //2.
129 fRowT0=0;
130 for (;(fNnodes+1)>(1<<fRowT0);fRowT0++);
131 fRowT0-=1;
132 // 2 = 2**0 + 1
133 // 3 = 2**1 + 1
134 // 4 = 2**1 + 2
135 // 5 = 2**2 + 1
136 // 6 = 2**2 + 2
137 // 7 = 2**2 + 3
138 // 8 = 2**2 + 4
139
140 //3.
141 // allocate space for boundaries
142 fRange = new Value[2*fNDim];
143 fIndPoints= new Index[fNpoints];
144 for (Index i=0; i<fNpoints; i++) fIndPoints[i] = i;
117c62ab 145 fAxis = new UChar_t[fNnodes];
146 fValue = new Value[fNnodes];
f2040a8f 147 //
148 fCrossNode = (1<<(fRowT0+1))-1;
149 if (fCrossNode<fNnodes) fCrossNode = 2*fCrossNode+1;
150 //
151 // fOffset = (((fNnodes+1)-(1<<fRowT0)))*2;
152 Int_t over = (fNnodes+1)-(1<<fRowT0);
153 Int_t filled = ((1<<fRowT0)-over)*fBucketSize;
154 fOffset = fNpoints-filled;
155 //
156// printf("Row0 %d\n", fRowT0);
157// printf("CrossNode %d\n", fCrossNode);
158// printf("Offset %d\n", fOffset);
159 //
160 //
161 //4.
162 // stack for non recursive build - size 128 bytes enough
163 Int_t rowStack[128];
164 Int_t nodeStack[128];
165 Int_t npointStack[128];
166 Int_t posStack[128];
167 Int_t currentIndex = 0;
168 Int_t iter =0;
169 rowStack[0] = 0;
170 nodeStack[0] = 0;
171 npointStack[0] = fNpoints;
172 posStack[0] = 0;
173 //
174 Int_t nbucketsall =0;
175 while (currentIndex>=0){
176 iter++;
177 //
178 Int_t npoints = npointStack[currentIndex];
179 if (npoints<=fBucketSize) {
180 //printf("terminal node : index %d iter %d\n", currentIndex, iter);
181 currentIndex--;
182 nbucketsall++;
183 continue; // terminal node
184 }
185 Int_t crow = rowStack[currentIndex];
186 Int_t cpos = posStack[currentIndex];
187 Int_t cnode = nodeStack[currentIndex];
f2040a8f 188 //printf("currentIndex %d npoints %d node %d\n", currentIndex, npoints, cnode);
189 //
190 // divide points
191 Int_t nbuckets0 = npoints/fBucketSize; //current number of buckets
192 if (npoints%fBucketSize) nbuckets0++; //
193 Int_t restRows = fRowT0-rowStack[currentIndex]; // rest of fully occupied node row
194 if (restRows<0) restRows =0;
195 for (;nbuckets0>(2<<restRows); restRows++);
196 Int_t nfull = 1<<restRows;
197 Int_t nrest = nbuckets0-nfull;
198 Int_t nleft =0, nright =0;
199 //
200 if (nrest>(nfull/2)){
201 nleft = nfull*fBucketSize;
202 nright = npoints-nleft;
203 }else{
204 nright = nfull*fBucketSize/2;
205 nleft = npoints-nright;
206 }
207
208 //
209 //find the axis with biggest spread
210 Value maxspread=0;
211 Value tempspread, min, max;
212 Index axspread=0;
213 Value *array;
214 for (Int_t idim=0; idim<fNDim; idim++){
215 array = fData[idim];
216 Spread(npoints, array, fIndPoints+cpos, min, max);
217 tempspread = max - min;
218 if (maxspread < tempspread) {
219 maxspread=tempspread;
220 axspread = idim;
221 }
117c62ab 222 if(cnode) continue;
223 //printf("set %d %6.3f %6.3f\n", idim, min, max);
224 fRange[2*idim] = min; fRange[2*idim+1] = max;
f2040a8f 225 }
226 array = fData[axspread];
227 KOrdStat(npoints, array, nleft, fIndPoints+cpos);
117c62ab 228 fAxis[cnode] = axspread;
229 fValue[cnode] = array[fIndPoints[cpos+nleft]];
f2040a8f 230 //printf("Set node %d : ax %d val %f\n", cnode, node->fAxis, node->fValue);
231 //
232 //
233 npointStack[currentIndex] = nleft;
234 rowStack[currentIndex] = crow+1;
235 posStack[currentIndex] = cpos;
236 nodeStack[currentIndex] = cnode*2+1;
237 currentIndex++;
238 npointStack[currentIndex] = nright;
239 rowStack[currentIndex] = crow+1;
240 posStack[currentIndex] = cpos+nleft;
241 nodeStack[currentIndex] = (cnode*2)+2;
242 //
243 if (0){
244 // consistency check
245 Info("Build()", Form("points %d left %d right %d", npoints, nleft, nright));
246 if (nleft<nright) Warning("Build", "Problem Left-Right");
247 if (nleft<0 || nright<0) Warning("Build()", "Problem Negative number");
248 }
249 }
f2040a8f 250
251 //printf("NBuckets\t%d\n", nbucketsall);
252 //fData = 0;
253}
254
255
256//_________________________________________________________________
257template <typename Index, typename Value>
117c62ab 258Bool_t TKDTree<Index, Value>::FindNearestNeighbors(const Value *point, const Int_t k, Index *&in, Value *&d)
f2040a8f 259{
260// Find "k" nearest neighbors to "point".
261//
262// Return true in case of success and false in case of failure.
263// The indexes of the nearest k points are stored in the array "in" in
264// increasing order of the distance to "point" and the maxim distance
265// in "d".
266//
117c62ab 267// The arrays "in" and "d" are managed internally by the TKDTree.
f2040a8f 268
269 Index inode = FindNode(point);
270 if(inode < fNnodes){
271 Error("FindNearestNeighbors()", "Point outside data range.");
272 return kFALSE;
273 }
274
316a7f5a 275 UInt_t debug = 0;
276 Int_t nCalculations = 0;
277
f2040a8f 278 // allocate working memory
316a7f5a 279 if(!fDistBuffer){
280 fDistBuffer = new Value[fBucketSize];
281 fIndBuffer = new Index[fBucketSize];
282 }
283 if(fkNNdim < k){
117c62ab 284 if(debug>=1){
316a7f5a 285 if(fkNN) printf("Reallocate memory %d -> %d \n", fkNNdim, 2*k);
286 else printf("Allocate %d memory\n", 2*k);
117c62ab 287 }
316a7f5a 288 fkNNdim = 2*k;
289 if(fkNN){
290 delete [] fkNN;
291 delete [] fkNNdist;
292 }
293 fkNN = new Index[fkNNdim];
294 fkNNdist = new Value[fkNNdim];
295 }
296 memset(fkNN, -1, k*sizeof(Index));
297 for(int i=0; i<k; i++) fkNNdist[i] = 9999.;
298 Index itmp, jtmp; Value ftmp, gtmp;
299
f2040a8f 300 // calculate number of boundaries with the data domain.
f2040a8f 301 if(!fBoundaries) MakeBoundaries();
302 Value *bounds = &fBoundaries[inode*2*fNDim];
117c62ab 303 Index nBounds = 0;
f2040a8f 304 for(int idim=0; idim<fNDim; idim++){
117c62ab 305 if((bounds[2*idim] - fRange[2*idim]) < 1.E-10) nBounds++;
306 if((bounds[2*idim+1] - fRange[2*idim+1]) < 1.E-10) nBounds++;
f2040a8f 307 }
316a7f5a 308 if(debug>=1) printf("Calculate boundaries [nBounds = %d]\n", nBounds);
117c62ab 309
310
f2040a8f 311 // traverse tree
117c62ab 312 UChar_t ax;
313 Value val, dif;
a3408ed3 314 Int_t nAllNodes = fNnodes + fNpoints/fBucketSize + ((fNpoints%fBucketSize)?1:0);
f2040a8f 315 Int_t nodeStack[128], nodeIn[128];
316 Index currentIndex = 0;
317 nodeStack[0] = inode;
318 nodeIn[0] = inode;
319 while(currentIndex>=0){
320 Int_t tnode = nodeStack[currentIndex];
321 Int_t entry = nodeIn[currentIndex];
322 currentIndex--;
316a7f5a 323 if(debug>=1) printf("Analyse tnode %d entry %d\n", tnode, entry);
324
f2040a8f 325 // check if node is still eligible
117c62ab 326 Int_t inode = (tnode-1)>>1; //tnode/2 + (tnode%2) - 1;
327 ax = fAxis[inode];
328 val = fValue[inode];
329 dif = point[ax] - val;
330 if((TMath::Abs(dif) > fkNNdist[k-1]) &&
331 ((dif > 0. && (tnode&1)) || (dif < 0. && !(tnode&1)))) {
332 if(debug>=1) printf("\tremove %d\n", (tnode-1)>>1);
f2040a8f 333
334 // mark bound
335 nBounds++;
336 // end all recursions
337 if(nBounds==2 * fNDim) break;
338 continue;
339 }
340
341 if(IsTerminal(tnode)){
316a7f5a 342 if(debug>=2) printf("\tProcess terminal node %d\n", tnode);
f2040a8f 343 // Link data to terminal node
344 Int_t offset = (tnode >= fCrossNode) ? (tnode-fCrossNode)*fBucketSize : fOffset+(tnode-fNnodes)*fBucketSize;
345 Index *indexPoints = &fIndPoints[offset];
316a7f5a 346 Int_t nbs = (tnode == 2*fNnodes) ? fNpoints%fBucketSize : fBucketSize;
347 nbs = nbs ? nbs : fBucketSize;
348 nCalculations += nbs;
349
350 Int_t npoints = 0;
351 for(int idx=0; idx<nbs; idx++){
f2040a8f 352 // calculate distance in the L1 metric
316a7f5a 353 fDistBuffer[npoints] = 0.;
354 for(int idim=0; idim<fNDim; idim++) fDistBuffer[npoints] += TMath::Abs(point[idim] - fData[idim][indexPoints[idx]]);
355 // register candidate neighbor
356 if(fDistBuffer[npoints] < fkNNdist[k-1]){
357 fIndBuffer[npoints] = indexPoints[idx];
358 npoints++;
359 }
f2040a8f 360 }
316a7f5a 361 for(int ibrowse=0; ibrowse<npoints; ibrowse++){
362 if(fDistBuffer[ibrowse] >= fkNNdist[k-1]) continue;
363 //insert neighbor
364 int iNN=0;
365 while(fDistBuffer[ibrowse] >= fkNNdist[iNN]) iNN++;
366 if(debug>=2) printf("\t\tinsert data %d @ %d distance %f\n", fIndBuffer[ibrowse], iNN, fDistBuffer[ibrowse]);
367
368 itmp = fkNN[iNN]; ftmp = fkNNdist[iNN];
369 fkNN[iNN] = fIndBuffer[ibrowse];
370 fkNNdist[iNN] = fDistBuffer[ibrowse];
371 for(int ireplace=iNN+1; ireplace<k; ireplace++){
372 jtmp = fkNN[ireplace]; gtmp = fkNNdist[ireplace];
373 fkNN[ireplace] = itmp; fkNNdist[ireplace] = ftmp;
374 itmp = jtmp; ftmp = gtmp;
375 if(ftmp == 9999.) break;
f2040a8f 376 }
316a7f5a 377 if(debug>=3){
378 for(int i=0; i<k; i++){
379 if(fkNNdist[i] == 9999.) break;
380 printf("\t\tfkNNdist[%d] = %f\n", i, fkNNdist[i]);
381 }
382 }
f2040a8f 383 }
384 }
385
386 // register parent
117c62ab 387 Int_t pn = (tnode-1)>>1;//tnode/2 + (tnode%2) - 1;
388 if(pn >= 0 && entry != pn){
f2040a8f 389 // check if parent node is eligible at all
117c62ab 390 if(TMath::Abs(point[fAxis[pn]] - fValue[pn]) > fkNNdist[k-1]){
f2040a8f 391 // mark bound
392 nBounds++;
393 // end all recursions
394 if(nBounds==2 * fNDim) break;
395 }
396 currentIndex++;
117c62ab 397 nodeStack[currentIndex]=pn;
f2040a8f 398 nodeIn[currentIndex]=tnode;
316a7f5a 399 if(debug>=2) printf("\tregister %d\n", nodeStack[currentIndex]);
f2040a8f 400 }
117c62ab 401 if(IsTerminal(tnode)) continue;
316a7f5a 402
f2040a8f 403 // register children nodes
117c62ab 404 Int_t cn;
f2040a8f 405 Bool_t kAllow[] = {kTRUE, kTRUE};
117c62ab 406 ax = fAxis[tnode];
407 val = fValue[tnode];
408 if(TMath::Abs(point[ax] - val) > fkNNdist[k-1]){
409 if(point[ax] > val) kAllow[0] = kFALSE;
f2040a8f 410 else kAllow[1] = kFALSE;
411 }
412 for(int ic=1; ic<=2; ic++){
413 if(!kAllow[ic-1]) continue;
117c62ab 414 cn = (tnode<<1)+ic;
415 if(cn < nAllNodes && entry != cn){
f2040a8f 416 currentIndex++;
117c62ab 417 nodeStack[currentIndex] = cn;
f2040a8f 418 nodeIn[currentIndex]=tnode;
316a7f5a 419 if(debug>=2) printf("\tregister %d\n", nodeStack[currentIndex]);
f2040a8f 420 }
316a7f5a 421 }
f2040a8f 422 }
423 // save results
424 in = fkNN;
316a7f5a 425 d = fkNNdist;
f2040a8f 426
117c62ab 427 return kTRUE;
f2040a8f 428}
429
430
431
432//_________________________________________________________________
433template <typename Index, typename Value>
434Index TKDTree<Index, Value>::FindNode(const Value * point){
435 //
436 // find the terminal node to which point belongs
117c62ab 437
f2040a8f 438 Index stackNode[128], inode;
439 Int_t currentIndex =0;
440 stackNode[0] = 0;
f2040a8f 441 while (currentIndex>=0){
442 inode = stackNode[currentIndex];
f2040a8f 443 if (IsTerminal(inode)) return inode;
444
117c62ab 445 currentIndex--;
446 if (point[fAxis[inode]]<=fValue[inode]){
f2040a8f 447 currentIndex++;
117c62ab 448 stackNode[currentIndex]=(inode<<1)+1;
f2040a8f 449 }
117c62ab 450 if (point[fAxis[inode]]>=fValue[inode]){
f2040a8f 451 currentIndex++;
117c62ab 452 stackNode[currentIndex]=(inode+1)<<1;
f2040a8f 453 }
454 }
a9c20b1f 455
456 return -1;
f2040a8f 457}
458
459
a9c20b1f 460
f2040a8f 461//_________________________________________________________________
462template <typename Index, typename Value>
463void TKDTree<Index, Value>::FindPoint(Value * point, Index &index, Int_t &iter){
464 //
465 // find the index of point
466 // works only if we keep fData pointers
467
468 Int_t stackNode[128];
469 Int_t currentIndex =0;
470 stackNode[0] = 0;
471 iter =0;
472 //
473 while (currentIndex>=0){
474 iter++;
475 Int_t inode = stackNode[currentIndex];
476 currentIndex--;
477 if (IsTerminal(inode)){
478 // investigate terminal node
479 Int_t indexIP = (inode >= fCrossNode) ? (inode-fCrossNode)*fBucketSize : (inode-fNnodes)*fBucketSize+fOffset;
480 printf("terminal %d indexP %d\n", inode, indexIP);
481 for (Int_t ibucket=0;ibucket<fBucketSize;ibucket++){
482 Bool_t isOK = kTRUE;
483 indexIP+=ibucket;
484 printf("ibucket %d index %d\n", ibucket, indexIP);
485 if (indexIP>=fNpoints) continue;
486 Int_t index0 = fIndPoints[indexIP];
487 for (Int_t idim=0;idim<fNDim;idim++) if (fData[idim][index0]!=point[idim]) isOK = kFALSE;
488 if (isOK) index = index0;
489 }
490 continue;
491 }
492
117c62ab 493 if (point[fAxis[inode]]<=fValue[inode]){
f2040a8f 494 currentIndex++;
495 stackNode[currentIndex]=(inode*2)+1;
496 }
117c62ab 497 if (point[fAxis[inode]]>=fValue[inode]){
f2040a8f 498 currentIndex++;
499 stackNode[currentIndex]=(inode*2)+2;
500 }
501 }
502 //
503 // printf("Iter\t%d\n",iter);
504}
505
506//_________________________________________________________________
507template <typename Index, typename Value>
508void TKDTree<Index, Value>::FindInRangeA(Value * point, Value * delta, Index *res , Index &npoints, Index & iter, Int_t bnode)
509{
510 //
511 // Find all points in the range specified by (point +- range)
512 // res - Resulting indexes are stored in res array
513 // npoints - Number of selected indexes in range
514 // NOTICE:
515 // For some cases it is better to don't keep data - because of memory consumption
516 // If the data are not kept - only boundary conditions are investigated
517 // some of the data can be outside of the range
518 // What is guranteed in this mode: All of the points in the range are selected + some fraction of others (but close)
519
520 Index stackNode[128];
521 iter=0;
522 Index currentIndex = 0;
523 stackNode[currentIndex] = bnode;
524 while (currentIndex>=0){
525 iter++;
526 Int_t inode = stackNode[currentIndex];
527 //
528 currentIndex--;
529 if (!IsTerminal(inode)){
530 // not terminal
117c62ab 531 //TKDNode<Index, Value> * node = &(fNodes[inode]);
532 if (point[fAxis[inode]] - delta[fAxis[inode]] < fValue[inode]) {
f2040a8f 533 currentIndex++;
534 stackNode[currentIndex]= (inode*2)+1;
535 }
117c62ab 536 if (point[fAxis[inode]] + delta[fAxis[inode]] >= fValue[inode]){
f2040a8f 537 currentIndex++;
538 stackNode[currentIndex]= (inode*2)+2;
539 }
540 }else{
541 Int_t indexIP =
542 (inode >= fCrossNode) ? (inode-fCrossNode)*fBucketSize : (inode-fNnodes)*fBucketSize+fOffset;
543 for (Int_t ibucket=0;ibucket<fBucketSize;ibucket++){
544 if (indexIP+ibucket>=fNpoints) break;
545 res[npoints] = fIndPoints[indexIP+ibucket];
546 npoints++;
547 }
548 }
549 }
550 if (fData){
551 //
552 // compress rest if data still accesible
553 //
554 Index npoints2 = npoints;
555 npoints=0;
556 for (Index i=0; i<npoints2;i++){
557 Bool_t isOK = kTRUE;
558 for (Index idim = 0; idim< fNDim; idim++){
559 if (TMath::Abs(fData[idim][res[i]]- point[idim])>delta[idim])
560 isOK = kFALSE;
561 }
562 if (isOK){
563 res[npoints] = res[i];
564 npoints++;
565 }
566 }
567 }
568}
569
570
571//_________________________________________________________________
572template <typename Index, typename Value>
573void TKDTree<Index, Value>::FindInRangeB(Value * point, Value * delta, Index *res , Index &npoints,Index & iter, Int_t bnode)
574{
575 //
576 Long64_t goldStatus = (1<<(2*fNDim))-1; // gold status
577 Index stackNode[128];
578 Long64_t stackStatus[128];
579 iter=0;
580 Index currentIndex = 0;
581 stackNode[currentIndex] = bnode;
582 stackStatus[currentIndex] = 0;
583 while (currentIndex>=0){
584 Int_t inode = stackNode[currentIndex];
585 Long64_t status = stackStatus[currentIndex];
586 currentIndex--;
587 iter++;
588 if (IsTerminal(inode)){
589 Int_t indexIP =
590 (inode >= fCrossNode) ? (inode-fCrossNode)*fBucketSize : (inode-fNnodes)*fBucketSize+fOffset;
591 for (Int_t ibucket=0;ibucket<fBucketSize;ibucket++){
592 if (indexIP+ibucket>=fNpoints) break;
593 res[npoints] = fIndPoints[indexIP+ibucket];
594 npoints++;
595 }
596 continue;
597 }
598 // not terminal
599 if (status == goldStatus){
600 Int_t ileft = inode;
601 Int_t iright = inode;
602 for (;ileft<fNnodes; ileft = (ileft<<1)+1);
603 for (;iright<fNnodes; iright = (iright<<1)+2);
604 Int_t indexL =
605 (ileft >= fCrossNode) ? (ileft-fCrossNode)*fBucketSize : (ileft-fNnodes)*fBucketSize+fOffset;
606 Int_t indexR =
607 (iright >= fCrossNode) ? (iright-fCrossNode)*fBucketSize : (iright-fNnodes)*fBucketSize+fOffset;
608 if (indexL<=indexR){
609 Int_t endpoint = indexR+fBucketSize;
610 if (endpoint>fNpoints) endpoint=fNpoints;
611 for (Int_t ipoint=indexL;ipoint<endpoint;ipoint++){
612 res[npoints] = fIndPoints[ipoint];
613 npoints++;
614 }
615 continue;
616 }
617 }
618 //
117c62ab 619 // TKDNode<Index, Value> * node = &(fNodes[inode]);
620 if (point[fAxis[inode]] - delta[fAxis[inode]] < fValue[inode]) {
f2040a8f 621 currentIndex++;
622 stackNode[currentIndex]= (inode<<1)+1;
117c62ab 623 if (point[fAxis[inode]] + delta[fAxis[inode]] > fValue[inode])
624 stackStatus[currentIndex]= status | (1<<(2*fAxis[inode]));
f2040a8f 625 }
117c62ab 626 if (point[fAxis[inode]] + delta[fAxis[inode]] >= fValue[inode]){
f2040a8f 627 currentIndex++;
628 stackNode[currentIndex]= (inode<<1)+2;
117c62ab 629 if (point[fAxis[inode]] - delta[fAxis[inode]]<fValue[inode])
630 stackStatus[currentIndex]= status | (1<<(2*fAxis[inode]+1));
f2040a8f 631 }
632 }
633 if (fData){
634 // compress rest
635 Index npoints2 = npoints;
636 npoints=0;
637 for (Index i=0; i<npoints2;i++){
638 Bool_t isOK = kTRUE;
639 for (Index idim = 0; idim< fNDim; idim++){
640 if (TMath::Abs(fData[idim][res[i]]- point[idim])>delta[idim])
641 isOK = kFALSE;
642 }
643 if (isOK){
644 res[npoints] = res[i];
645 npoints++;
646 }
647 }
648 }
649}
650
651
652
653//_________________________________________________________________
654template <typename Index, typename Value>
655void TKDTree<Index, Value>::SetData(Index npoints, Index ndim, UInt_t bsize, Value **data)
656{
657// TO DO
658//
659// Check reconstruction/reallocation of memory of data. Maybe it is not
660// necessary to delete and realocate space but only to use the same space
661
662 Clear();
663
664 //Columnwise!!!!
665 fData = data;
666 fNpoints = npoints;
667 fNDim = ndim;
668 fBucketSize = bsize;
669
670 Build();
671}
672
673
674
675//_________________________________________________________________
676template <typename Index, typename Value>
677void TKDTree<Index, Value>::Spread(Index ntotal, Value *a, Index *index, Value &min, Value &max)
678{
679 //Value min, max;
680 Index i;
681 min = a[index[0]];
682 max = a[index[0]];
683 for (i=0; i<ntotal; i++){
684 if (a[index[i]]<min) min = a[index[i]];
685 if (a[index[i]]>max) max = a[index[i]];
686 }
687}
688
689
690//_________________________________________________________________
691template <typename Index, typename Value>
692Value TKDTree<Index, Value>::KOrdStat(Index ntotal, Value *a, Index k, Index *index)
693{
694 //
695 //copy of the TMath::KOrdStat because I need an Index work array
696
697 Index i, ir, j, l, mid;
698 Index arr;
699 Index temp;
700
701 Index rk = k;
702 l=0;
703 ir = ntotal-1;
704 for(;;) {
705 if (ir<=l+1) { //active partition contains 1 or 2 elements
706 if (ir == l+1 && a[index[ir]]<a[index[l]])
707 {temp = index[l]; index[l]=index[ir]; index[ir]=temp;}
708 Value tmp = a[index[rk]];
709 return tmp;
710 } else {
711 mid = (l+ir) >> 1; //choose median of left, center and right
712 {temp = index[mid]; index[mid]=index[l+1]; index[l+1]=temp;}//elements as partitioning element arr.
713 if (a[index[l]]>a[index[ir]]) //also rearrange so that a[l]<=a[l+1]
714 {temp = index[l]; index[l]=index[ir]; index[ir]=temp;}
715
716 if (a[index[l+1]]>a[index[ir]])
717 {temp=index[l+1]; index[l+1]=index[ir]; index[ir]=temp;}
718
719 if (a[index[l]]>a[index[l+1]])
720 {temp = index[l]; index[l]=index[l+1]; index[l+1]=temp;}
721
722 i=l+1; //initialize pointers for partitioning
723 j=ir;
724 arr = index[l+1];
725 for (;;) {
726 do i++; while (a[index[i]]<a[arr]);
727 do j--; while (a[index[j]]>a[arr]);
728 if (j<i) break; //pointers crossed, partitioning complete
729 {temp=index[i]; index[i]=index[j]; index[j]=temp;}
730 }
731 index[l+1]=index[j];
732 index[j]=arr;
733 if (j>=rk) ir = j-1; //keep active the partition that
734 if (j<=rk) l=i; //contains the k_th element
735 }
736 }
737}
738
739//_________________________________________________________________
740template <typename Index, typename Value>
741void TKDTree<Index, Value>::MakeBoundaries(Value *range)
742{
743// Build boundaries for each node.
744
745
117c62ab 746 if(range) memcpy(fRange, range, fNDimm*sizeof(Value));
f2040a8f 747 // total number of nodes including terminal nodes
a3408ed3 748 Int_t totNodes = fNnodes + fNpoints/fBucketSize + ((fNpoints%fBucketSize)?1:0);
117c62ab 749 fBoundaries = new Value[totNodes*fNDimm];
750 //Info("MakeBoundaries(Value*)", Form("Allocate boundaries for %d nodes", totNodes));
751
316a7f5a 752
f2040a8f 753 // loop
117c62ab 754 Value *tbounds = 0x0, *cbounds = 0x0;
755 Int_t cn;
f2040a8f 756 for(int inode=fNnodes-1; inode>=0; inode--){
117c62ab 757 tbounds = &fBoundaries[inode*fNDimm];
758 memcpy(tbounds, fRange, fNDimm*sizeof(Value));
f2040a8f 759
760 // check left child node
117c62ab 761 cn = (inode<<1)+1;
762 if(IsTerminal(cn)) CookBoundaries(inode, kTRUE);
763 cbounds = &fBoundaries[fNDimm*cn];
764 for(int idim=0; idim<fNDim; idim++) tbounds[idim<<1] = cbounds[idim<<1];
f2040a8f 765
766 // check right child node
117c62ab 767 cn = (inode+1)<<1;
768 if(IsTerminal(cn)) CookBoundaries(inode, kFALSE);
769 cbounds = &fBoundaries[fNDimm*cn];
770 for(int idim=0; idim<fNDim; idim++) tbounds[(idim<<1)+1] = cbounds[(idim<<1)+1];
f2040a8f 771 }
772}
773
774//_________________________________________________________________
775template <typename Index, typename Value>
117c62ab 776void TKDTree<Index, Value>::CookBoundaries(const Int_t node, Bool_t LEFT)
f2040a8f 777{
778 // define index of this terminal node
117c62ab 779 Int_t index = (node<<1) + (LEFT ? 1 : 2);
780 //Info("CookBoundaries()", Form("Node %d", index));
f2040a8f 781
782 // build and initialize boundaries for this node
117c62ab 783 Value *tbounds = &fBoundaries[fNDimm*index];
784 memcpy(tbounds, fRange, fNDimm*sizeof(Value));
f2040a8f 785 Bool_t flag[256]; // cope with up to 128 dimensions
117c62ab 786 memset(flag, kFALSE, fNDimm);
f2040a8f 787 Int_t nvals = 0;
788
789 // recurse parent nodes
117c62ab 790 Int_t pn = node;
791 while(pn >= 0 && nvals < fNDimm){
f2040a8f 792 if(LEFT){
117c62ab 793 index = (fAxis[pn]<<1)+1;
794 if(!flag[index]) {
795 tbounds[index] = fValue[pn];
796 flag[index] = kTRUE;
f2040a8f 797 nvals++;
798 }
799 } else {
117c62ab 800 index = fAxis[pn]<<1;
801 if(!flag[index]) {
802 tbounds[index] = fValue[pn];
803 flag[index] = kTRUE;
f2040a8f 804 nvals++;
805 }
806 }
117c62ab 807 LEFT = pn&1;
808 pn = (pn - 1)>>1;
f2040a8f 809 }
810}
811
f2040a8f 812template class TKDTree<Int_t, Float_t>;
813template class TKDTree<Int_t, Double_t>;
814