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