775c6022fe400ec09977c9098d1e72ecf712870d
[u/mrichter/AliRoot.git] / STAT / TKDTree.cxx
1 #include "TKDTree.h"
2
3 #include "TString.h"
4 #include <string.h>
5
6 #ifndef R__ALPHA
7 templateClassImp(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 ///////////////////////////////////////////////////////////////////////
32
33
34 //_________________________________________________________________
35 template <typename  Index, typename Value>
36 TKDTree<Index, Value>::TKDTree() :
37         TObject()
38         ,kDataOwner(kFALSE)
39         ,fNnodes(0)
40         ,fNDim(0)
41         ,fNDimm(0)
42         ,fNpoints(0)
43         ,fBucketSize(0)
44         ,fAxis(0x0)
45         ,fValue(0x0)
46         ,fRange(0x0)
47         ,fData(0x0)
48         ,fBoundaries(0x0)
49         ,fkNNdim(0)
50         ,fkNN(0x0)
51         ,fkNNdist(0x0)
52         ,fDistBuffer(0x0)
53         ,fIndBuffer(0x0)
54         ,fIndPoints(0x0)
55         ,fRowT0(0)
56         ,fCrossNode(0)
57         ,fOffset(0)
58 {
59 // Default constructor. Nothing is built
60 }
61
62
63 //_________________________________________________________________
64 template <typename  Index, typename Value>
65 TKDTree<Index, Value>::TKDTree(Index npoints, Index ndim, UInt_t bsize, Value **data) :
66         TObject()
67         ,kDataOwner(kFALSE)
68         ,fNnodes(0)
69         ,fNDim(ndim)
70         ,fNDimm(2*ndim)
71         ,fNpoints(npoints)
72         ,fBucketSize(bsize)
73         ,fAxis(0x0)
74         ,fValue(0x0)
75         ,fRange(0x0)
76         ,fData(data) //Columnwise!!!!!
77         ,fBoundaries(0x0)
78         ,fkNNdim(0)
79         ,fkNN(0x0)
80         ,fkNNdist(0x0)
81         ,fDistBuffer(0x0)
82         ,fIndBuffer(0x0)
83         ,fIndPoints(0x0)
84         ,fRowT0(0)
85         ,fCrossNode(0)
86         ,fOffset(0)
87 {
88 // Allocate data by hand. See TKDTree(TTree*, const Char_t*) for an automatic method.
89
90         Build();
91 }
92
93 //_________________________________________________________________
94 template <typename  Index, typename Value>
95 TKDTree<Index, Value>::~TKDTree()
96 {
97         if (fIndBuffer) delete [] fIndBuffer;
98         if (fDistBuffer) delete [] fDistBuffer;
99         if (fkNNdist) delete [] fkNNdist;
100         if (fkNN) delete [] fkNN;
101         if (fAxis) delete [] fAxis;
102         if (fValue) delete [] fValue;
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 //_________________________________________________________________
114 template <typename  Index, typename Value>
115 void 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;
145         fAxis  = new UChar_t[fNnodes];
146         fValue = new Value[fNnodes];
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];               
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                         }
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;
225                 }
226                 array = fData[axspread];
227                 KOrdStat(npoints, array, nleft, fIndPoints+cpos);
228                 fAxis[cnode]  = axspread;
229                 fValue[cnode] = array[fIndPoints[cpos+nleft]];
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         }
250         
251         //printf("NBuckets\t%d\n", nbucketsall);
252         //fData = 0;
253 }
254
255
256 //_________________________________________________________________
257 template <typename  Index, typename Value>
258 Bool_t  TKDTree<Index, Value>::FindNearestNeighbors(const Value *point, const Int_t k, Index *&in, Value *&d)
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 //
267 // The arrays "in" and "d" are managed internally by the TKDTree.
268
269         Index inode = FindNode(point);
270         if(inode < fNnodes){
271                 Error("FindNearestNeighbors()", "Point outside data range.");
272                 return kFALSE;
273         }
274
275         UInt_t debug = 0;
276         Int_t nCalculations = 0;
277         
278         // allocate working memory
279         if(!fDistBuffer){
280                 fDistBuffer = new Value[fBucketSize];
281                 fIndBuffer  = new Index[fBucketSize];
282         }
283         if(fkNNdim < k){
284                 if(debug>=1){
285                         if(fkNN) printf("Reallocate memory %d -> %d \n", fkNNdim, 2*k);
286                         else printf("Allocate %d memory\n", 2*k);
287                 }
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         
300         // calculate number of boundaries with the data domain. 
301         if(!fBoundaries) MakeBoundaries();
302         Value *bounds = &fBoundaries[inode*2*fNDim];
303         Index nBounds = 0;
304         for(int idim=0; idim<fNDim; idim++){
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++;
307         }
308         if(debug>=1) printf("Calculate boundaries [nBounds = %d]\n", nBounds);
309
310                 
311         // traverse tree
312         UChar_t ax;
313         Value   val, dif;
314         Int_t nAllNodes = fNnodes + fNpoints/fBucketSize + ((fNpoints%fBucketSize)?1:0);
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--;
323                 if(debug>=1) printf("Analyse tnode %d entry %d\n", tnode, entry);
324                 
325                 // check if node is still eligible
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);
333
334                         // mark bound
335                         nBounds++;
336                         // end all recursions
337                         if(nBounds==2 * fNDim) break;
338                         continue;
339                 }
340                 
341                 if(IsTerminal(tnode)){
342                         if(debug>=2) printf("\tProcess terminal node %d\n", tnode);
343                         // Link data to terminal node
344                         Int_t offset = (tnode >= fCrossNode) ? (tnode-fCrossNode)*fBucketSize : fOffset+(tnode-fNnodes)*fBucketSize;
345                         Index *indexPoints = &fIndPoints[offset];
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++){
352                                 // calculate distance in the L1 metric
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                                 }
360                         }
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;
376                                 }
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                                 }                               
383                         }
384                 }
385                 
386                 // register parent
387                 Int_t pn = (tnode-1)>>1;//tnode/2 + (tnode%2) - 1;
388                 if(pn >= 0 && entry != pn){
389                         // check if parent node is eligible at all
390                         if(TMath::Abs(point[fAxis[pn]] - fValue[pn]) > fkNNdist[k-1]){
391                                 // mark bound
392                                 nBounds++;
393                                 // end all recursions
394                                 if(nBounds==2 * fNDim) break;
395                         }
396                         currentIndex++;
397                         nodeStack[currentIndex]=pn;
398                         nodeIn[currentIndex]=tnode;
399                         if(debug>=2) printf("\tregister %d\n", nodeStack[currentIndex]);
400                 }
401                 if(IsTerminal(tnode)) continue;
402                 
403                 // register children nodes
404                 Int_t cn;
405                 Bool_t kAllow[] = {kTRUE, kTRUE};
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;
410                         else kAllow[1] = kFALSE;
411                 }
412                 for(int ic=1; ic<=2; ic++){
413                         if(!kAllow[ic-1]) continue;
414                         cn = (tnode<<1)+ic;
415                         if(cn < nAllNodes && entry != cn){
416                                 currentIndex++;
417                                 nodeStack[currentIndex] = cn;
418                                 nodeIn[currentIndex]=tnode;
419                                 if(debug>=2) printf("\tregister %d\n", nodeStack[currentIndex]);
420                         }
421                 }               
422         }
423         // save results
424         in = fkNN;
425         d  = fkNNdist;
426         
427         return kTRUE;
428 }
429
430
431
432 //_________________________________________________________________
433 template <typename  Index, typename Value>
434 Index TKDTree<Index, Value>::FindNode(const Value * point){
435   //
436   // find the terminal node to which point belongs
437
438         Index stackNode[128], inode;
439         Int_t currentIndex =0;
440         stackNode[0] = 0;
441         while (currentIndex>=0){
442                 inode    = stackNode[currentIndex];
443                 if (IsTerminal(inode)) return inode;
444                 
445                 currentIndex--;
446                 if (point[fAxis[inode]]<=fValue[inode]){
447                         currentIndex++;
448                         stackNode[currentIndex]=(inode<<1)+1;
449                 }
450                 if (point[fAxis[inode]]>=fValue[inode]){
451                         currentIndex++;
452                         stackNode[currentIndex]=(inode+1)<<1;
453                 }
454         }
455         
456         return -1;
457 }
458
459
460
461 //_________________________________________________________________
462 template <typename  Index, typename Value>
463 void 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                 
493                 if (point[fAxis[inode]]<=fValue[inode]){
494                         currentIndex++;
495                         stackNode[currentIndex]=(inode*2)+1;
496                 }
497                 if (point[fAxis[inode]]>=fValue[inode]){
498                         currentIndex++;
499                         stackNode[currentIndex]=(inode*2)+2;
500                 }
501         }
502   //
503   //  printf("Iter\t%d\n",iter);
504 }
505
506 //_________________________________________________________________
507 template <typename  Index, typename Value>
508 void 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
531       //TKDNode<Index, Value> * node = &(fNodes[inode]);
532       if (point[fAxis[inode]] - delta[fAxis[inode]] <  fValue[inode]) {
533         currentIndex++; 
534         stackNode[currentIndex]= (inode*2)+1;
535       }
536       if (point[fAxis[inode]] + delta[fAxis[inode]] >= fValue[inode]){
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 //_________________________________________________________________
572 template <typename  Index, typename Value>
573 void 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     //
619     // TKDNode<Index, Value> * node = &(fNodes[inode]);
620     if (point[fAxis[inode]] - delta[fAxis[inode]] <  fValue[inode]) {
621       currentIndex++; 
622       stackNode[currentIndex]= (inode<<1)+1;
623       if (point[fAxis[inode]] + delta[fAxis[inode]] > fValue[inode])
624         stackStatus[currentIndex]= status | (1<<(2*fAxis[inode]));
625     }
626     if (point[fAxis[inode]] + delta[fAxis[inode]] >= fValue[inode]){
627       currentIndex++; 
628       stackNode[currentIndex]= (inode<<1)+2;
629       if (point[fAxis[inode]] - delta[fAxis[inode]]<fValue[inode])
630         stackStatus[currentIndex]= status | (1<<(2*fAxis[inode]+1));
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 //_________________________________________________________________
654 template <typename  Index, typename Value>
655 void 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 //_________________________________________________________________
676 template <typename  Index, typename Value>
677 void 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 //_________________________________________________________________
691 template <typename  Index, typename Value>
692 Value 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 //_________________________________________________________________
740 template <typename Index, typename Value>
741 void TKDTree<Index, Value>::MakeBoundaries(Value *range)
742 {
743 // Build boundaries for each node.
744
745
746         if(range) memcpy(fRange, range, fNDimm*sizeof(Value));
747         // total number of nodes including terminal nodes
748         Int_t totNodes = fNnodes + fNpoints/fBucketSize + ((fNpoints%fBucketSize)?1:0);
749         fBoundaries = new Value[totNodes*fNDimm];
750         //Info("MakeBoundaries(Value*)", Form("Allocate boundaries for %d nodes", totNodes));
751
752         
753         // loop
754         Value *tbounds = 0x0, *cbounds = 0x0;
755         Int_t cn;
756         for(int inode=fNnodes-1; inode>=0; inode--){
757                 tbounds = &fBoundaries[inode*fNDimm];
758                 memcpy(tbounds, fRange, fNDimm*sizeof(Value));
759                                 
760                 // check left child node
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];
765                 
766                 // check right child node
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];
771         }
772 }
773
774 //_________________________________________________________________
775 template <typename Index, typename Value>
776 void TKDTree<Index, Value>::CookBoundaries(const Int_t node, Bool_t LEFT)
777 {
778         // define index of this terminal node
779         Int_t index = (node<<1) + (LEFT ? 1 : 2);
780         //Info("CookBoundaries()", Form("Node %d", index));
781         
782         // build and initialize boundaries for this node
783         Value *tbounds = &fBoundaries[fNDimm*index];
784         memcpy(tbounds, fRange, fNDimm*sizeof(Value));
785         Bool_t flag[256];  // cope with up to 128 dimensions 
786         memset(flag, kFALSE, fNDimm);
787         Int_t nvals = 0;
788                 
789         // recurse parent nodes
790         Int_t pn = node;
791         while(pn >= 0 && nvals < fNDimm){
792                 if(LEFT){
793                         index = (fAxis[pn]<<1)+1;
794                         if(!flag[index]) {
795                                 tbounds[index] = fValue[pn];
796                                 flag[index] = kTRUE;
797                                 nvals++;
798                         }
799                 } else {
800                         index = fAxis[pn]<<1;
801                         if(!flag[index]) {
802                                 tbounds[index] = fValue[pn];
803                                 flag[index] = kTRUE;
804                                 nvals++;
805                         }
806                 }
807                 LEFT = pn&1;
808                 pn =  (pn - 1)>>1;
809         }
810 }
811
812 template class TKDTree<Int_t, Float_t>;
813 template class TKDTree<Int_t, Double_t>;
814