Logo ROOT   6.08/07
Reference Guide
TKDTreeBinning.cxx
Go to the documentation of this file.
1 // @(#)root/mathcore:$Id$
2 // Authors: B. Rabacal 11/2010
3 
4 /**********************************************************************
5  * *
6  * Copyright (c) 2010 , LCG ROOT MathLib Team *
7  * *
8  * *
9  **********************************************************************/
10 
11 // implementation file for class TKDTreeBinning
12 //
13 
14 #include <algorithm>
15 #include <limits>
16 #include <cmath>
17 
18 #include "TKDTreeBinning.h"
19 
20 #include "Fit/BinData.h"
21 #include "TRandom.h"
22 #include "TBuffer.h"
23 
25 
26 /**
27 \class TKDTreeBinning
28 <- TKDTreeBinning - A class providing multidimensional binning ->
29 
30 The class implements multidimensional binning by constructing a TKDTree inner structure from the
31 data which is used as the bins.
32 The bins are retrieved as two double*, one for the minimum bin edges,
33 the other as the maximum bin edges. For one dimension one of these is enough to correctly define the bins.
34 For the multidimensional case both minimum and maximum ones are necessary for the bins to be well defined.
35 The bin edges of d-dimensional data is a d-tet of the bin's thresholds. For example if d=3 the minimum bin
36 edges of bin b is of the form of the following array: {xbmin, ybmin, zbmin}.
37 You also have the possibility to sort the bins by their density.
38 
39 Details of usage can be found in `$ROOTSYS/tutorials/math/kdTreeBinning.C` and more information on
40 the embedded TKDTree documentation.
41 
42 @ingroup MathCore
43 
44 */
45 
47  // Boolean functor whose predicate depends on the bin's density. Used for ascending sort.
48  CompareAsc(const TKDTreeBinning* treebins) : bins(treebins) {}
49  Bool_t operator()(UInt_t bin1, UInt_t bin2) {
50  return bins->GetBinDensity(bin1) < bins->GetBinDensity(bin2);
51  }
52  const TKDTreeBinning* bins;
53 };
54 
56  // Boolean functor whose predicate depends on the bin's density. Used for descending sort.
57  CompareDesc(const TKDTreeBinning* treebins) : bins(treebins) {}
58  Bool_t operator()(UInt_t bin1, UInt_t bin2) {
59  return bins->GetBinDensity(bin1) > bins->GetBinDensity(bin2);
60  }
61  const TKDTreeBinning* bins;
62 };
63 
64 /// Class's constructor taking the size of the data points, dimension, a data array and the number
65 /// of bins (default = 100). It is reccomended to have the number of bins as an exact divider of
66 /// the data size.
67 /// The data array must be organized with a stride=1 for the points and = N (the dataSize) for the dimension.
68 ///
69 /// Thus data[] = x1,x2,x3,......xN, y1,y2,y3......yN, z1,z2,...........zN,....
70 ///
71 /// Note that the passed dataSize is not the size of the array but is the number of points (N)
72 /// The size of the array must be at least dataDim*dataSize
73 ///
74 TKDTreeBinning::TKDTreeBinning(UInt_t dataSize, UInt_t dataDim, Double_t* data, UInt_t nBins, bool adjustBinEdges)
75 : fData(0), fBinMinEdges(std::vector<Double_t>()), fBinMaxEdges(std::vector<Double_t>()), fDataBins((TKDTreeID*)0), fDim(dataDim),
76 fDataSize(dataSize), fDataThresholds(std::vector<std::pair<Double_t, Double_t> >(fDim, std::make_pair(0., 0.))),
78  if (adjustBinEdges) SetBit(kAdjustBinEdges);
79  if (data) {
80  SetData(data);
81  SetNBins(nBins);
82  } else {
83  if (fData.empty())
84  this->Warning("TKDTreeBinning", "Data is nil. Nothing is built.");
85  }
86 }
87 /// Class's constructor taking the size of the data points, dimension, a data vector and the number
88 /// of bins (default = 100). It is reccomended to have the number of bins as an exact divider of
89 /// the data size.
90 /// The data array must be organized with a stride=1 for the points and = N (the dataSize) for the dimension.
91 ///
92 /// Thus data[] = x1,x2,x3,......xN, y1,y2,y3......yN, z1,z2,...........zN,....
93 ///
94 /// Note that the passed data vector may contains a larger size, in case extra coordinates are associated but not used
95 /// in building the kdtree
96 /// The size of thedata vector must be at least dataDim*dataSize
97 ///
98 TKDTreeBinning::TKDTreeBinning(UInt_t dataSize, UInt_t dataDim, const std::vector<double> &data, UInt_t nBins, bool adjustBinEdges)
99 : fData(0), fBinMinEdges(std::vector<Double_t>()), fBinMaxEdges(std::vector<Double_t>()), fDataBins((TKDTreeID*)0), fDim(dataDim),
100 fDataSize(dataSize), fDataThresholds(std::vector<std::pair<Double_t, Double_t> >(fDim, std::make_pair(0., 0.))),
102  if (adjustBinEdges) SetBit(kAdjustBinEdges);
103  if (!data.empty()) {
104  SetData(data);
105  SetNBins(nBins);
106  } else {
107  if (fData.empty())
108  this->Warning("TKDTreeBinning", "Data is nil. Nothing is built.");
109  }
110 }
111 
112 /// Default constructor (for I/O)
114 // fData(nullptr),
116  fDim(0),
117  fDataSize(0),
119 {}
120 
121 /// Class's destructor
123  // if (fData) delete[] fData;
124  if (fDataBins) delete fDataBins;
125 }
126 
127 /// Sets binning inner structure
129  fNBins = bins;
130  if (fDim && fNBins && fDataSize) {
131  if (fDataSize / fNBins) {
132  Bool_t remainingData = fDataSize % fNBins;
133  if (remainingData) {
134  fNBins += 1;
135  this->Info("SetNBins", "Number of bins is not enough to hold the data. Extra bin added.");
136  }
137  fDataBins = new TKDTreeID(fDataSize, fDim, fDataSize / (fNBins - remainingData)); // TKDTree input is data size, data dimension and the content size of bins ("bucket" size)
138  SetTreeData();
139  fDataBins->Build();
140  SetBinsEdges();
141  SetBinsContent();
142  } else {
143  fDataBins = (TKDTreeID*)0;
144  this->Warning("SetNBins", "Number of bins is bigger than data size. Nothing is built.");
145  }
146  } else {
147  fDataBins = (TKDTreeID*)0;
148  if (!fDim)
149  this->Warning("SetNBins", "Data dimension is nil. Nothing is built.");
150  if (!fNBins)
151  this->Warning("SetNBins", "Number of bins is nil. Nothing is built.");
152  if (!fDataSize)
153  this->Warning("SetNBins", "Data size is nil. Nothing is built.");
154  }
155 }
156 
157 /// Sorts bins by their density
159  if (fDim == 1) {
160  // in one dim they are already sorted (no need to do anything)
161  return;
162  } else {
163  std::vector<UInt_t> indices(fNBins); // vector for indices (for the inverse transformation)
164  for (UInt_t i = 0; i < fNBins; ++i)
165  indices[i] = i;
166  if (sortAsc) {
167  std::sort(indices.begin(), indices.end(), CompareAsc(this));
169  } else {
170  std::sort(indices.begin(), indices.end(), CompareDesc(this));
172  }
173  std::vector<Double_t> binMinEdges(fNBins * fDim);
174  std::vector<Double_t> binMaxEdges(fNBins * fDim);
175  std::vector<UInt_t> binContent(fNBins ); // reajust also content (not needed bbut better in general!)
176  fIndices.resize(fNBins);
177  for (UInt_t i = 0; i < fNBins; ++i) {
178  for (UInt_t j = 0; j < fDim; ++j) {
179  binMinEdges[i * fDim + j] = fBinMinEdges[indices[i] * fDim + j];
180  binMaxEdges[i * fDim + j] = fBinMaxEdges[indices[i] * fDim + j];
181  }
182  binContent[i] = fBinsContent[indices[i]];
183  fIndices[indices[i]] = i;
184  }
185  fBinMinEdges.swap(binMinEdges);
186  fBinMaxEdges.swap(binMaxEdges);
187  fBinsContent.swap(binContent);
188 
189  // not needed anymore if readjusting bin content all the time
190  // re-adjust content of extra bins if exists
191  // since it is different than the others
192  // if ( fDataSize % fNBins != 0) {
193  // UInt_t k = 0;
194  // Bool_t found = kFALSE;
195  // while (!found) {
196  // if (indices[k] == fNBins - 1) {
197  // found = kTRUE;
198  // break;
199  // }
200  // ++k;
201  // }
202  // fBinsContent[fNBins - 1] = fDataBins->GetBucketSize();
203  // fBinsContent[k] = fDataSize % fNBins-1;
204  // }
205 
206  fIsSorted = kTRUE;
207  }
208 }
209 
211  // Sets the data and finds minimum and maximum by dimensional coordinate
212  fData.resize(fDim*fDataSize);
213  auto first = fData.begin();
214  for (UInt_t i = 0; i < fDim; ++i) {
215  for (UInt_t j = 0; j < fDataSize; ++j) {
216  fData[i*fDataSize+j] = data[i * fDataSize + j];
217  }
218  auto end = first+fDataSize;
219  fDataThresholds[i] = std::make_pair(*std::min_element(first, end), *std::max_element(first,end));
220  first = end;
221  }
222 }
223 void TKDTreeBinning::SetData(const std::vector<double>& data) {
224  // Sets the data and finds minimum and maximum by dimensional coordinate
225  fData = data;
226  auto first = fData.begin();
227  // find min/max
228  for (UInt_t i = 0; i < fDim; ++i) {
229  auto end = first+fDataSize;
230  fDataThresholds[i] = std::make_pair(*std::min_element(first, end), *std::max_element(first,end));
231  first = end;
232  }
233 }
234 
236  // Sets the data for constructing the kD-tree
237  for (UInt_t i = 0; i < fDim; ++i)
238  fDataBins->SetData(i, &fData[i*fDataSize]);
239 }
240 
242  // Sets the bins' content
243  fBinsContent.reserve(fNBins);
244  for (UInt_t i = 0; i < fNBins; ++i)
246  if ( fDataSize % fNBins != 0 )
247  fBinsContent[fNBins - 1] = fDataSize % (fNBins-1);
248 }
249 
251  // Sets the bins' edges
252  //Double_t* rawBinEdges = fDataBins->GetBoundaryExact(fDataBins->GetNNodes());
253  Double_t* rawBinEdges = fDataBins->GetBoundary(fDataBins->GetNNodes());
254  fCheckedBinEdges = std::vector<std::vector<std::pair<Bool_t, Bool_t> > >(fDim, std::vector<std::pair<Bool_t, Bool_t> >(fNBins, std::make_pair(kFALSE, kFALSE)));
255  fCommonBinEdges = std::vector<std::map<Double_t, std::vector<UInt_t> > >(fDim, std::map<Double_t, std::vector<UInt_t> >());
256  SetCommonBinEdges(rawBinEdges);
257  if (TestBit(kAdjustBinEdges) ) {
258  ReadjustMinBinEdges(rawBinEdges);
259  ReadjustMaxBinEdges(rawBinEdges);
260  }
261  SetBinMinMaxEdges(rawBinEdges);
262  fCommonBinEdges.clear();
263  fCheckedBinEdges.clear();
264 }
265 
267  // Sets the bins' minimum and maximum edges
268  fBinMinEdges.reserve(fNBins * fDim);
269  fBinMaxEdges.reserve(fNBins * fDim);
270  for (UInt_t i = 0; i < fNBins; ++i) {
271  for (UInt_t j = 0; j < fDim; ++j) {
272  fBinMinEdges.push_back(binEdges[(i * fDim + j) * 2]);
273  fBinMaxEdges.push_back(binEdges[(i * fDim + j) * 2 + 1]);
274  }
275  }
276 }
277 
279  // Sets indexing on the bin edges which have common boundaries
280  for (UInt_t i = 0; i < fDim; ++i) {
281  for (UInt_t j = 0; j < fNBins; ++j) {
282  Double_t binEdge = binEdges[(j * fDim + i) * 2];
283  if(fCommonBinEdges[i].find(binEdge) == fCommonBinEdges[i].end()) {
284  std::vector<UInt_t> commonBinEdges;
285  for (UInt_t k = 0; k < fNBins; ++k) {
286  UInt_t minBinEdgePos = (k * fDim + i) * 2;
287  if (std::fabs(binEdge - binEdges[minBinEdgePos]) < std::numeric_limits<Double_t>::epsilon())
288  commonBinEdges.push_back(minBinEdgePos);
289  UInt_t maxBinEdgePos = ++minBinEdgePos;
290  if (std::fabs(binEdge - binEdges[maxBinEdgePos]) < std::numeric_limits<Double_t>::epsilon())
291  commonBinEdges.push_back(maxBinEdgePos);
292  }
293  fCommonBinEdges[i][binEdge] = commonBinEdges;
294  }
295  }
296  }
297 }
298 
300  // Readjusts the bins' minimum edge by shifting it slightly lower
301  // to avoid overlapping with the data
302  for (UInt_t i = 0; i < fDim; ++i) {
303  for (UInt_t j = 0; j < fNBins; ++j) {
304  if (!fCheckedBinEdges[i][j].first) {
305  Double_t binEdge = binEdges[(j * fDim + i) * 2];
306  Double_t adjustedBinEdge = binEdge;
307  double eps = -10*std::numeric_limits<Double_t>::epsilon();
308  if (adjustedBinEdge != 0)
309  adjustedBinEdge *= (1. + eps);
310  else
311  adjustedBinEdge += eps;
312 
313  for (UInt_t k = 0; k < fCommonBinEdges[i][binEdge].size(); ++k) {
314  UInt_t binEdgePos = fCommonBinEdges[i][binEdge][k];
315  Bool_t isMinBinEdge = binEdgePos % 2 == 0;
316  UInt_t bin = isMinBinEdge ? (binEdgePos / 2 - i) / fDim : ((binEdgePos - 1) / 2 - i) / fDim;
317  binEdges[binEdgePos] = adjustedBinEdge;
318  if (isMinBinEdge)
319  fCheckedBinEdges[i][bin].first = kTRUE;
320  else
321  fCheckedBinEdges[i][bin].second = kTRUE;
322  }
323  }
324  }
325  }
326 }
327 
329  // Readjusts the bins' maximum edge
330  // and shift it sligtly higher
331  for (UInt_t i = 0; i < fDim; ++i) {
332  for (UInt_t j = 0; j < fNBins; ++j) {
333  if (!fCheckedBinEdges[i][j].second) {
334  Double_t& binEdge = binEdges[(j * fDim + i) * 2 + 1];
335  double eps = 10*std::numeric_limits<Double_t>::epsilon();
336  if (binEdge != 0)
337  binEdge *= (1. + eps);
338  else
339  binEdge += eps;
340 
341 
342  }
343  }
344  }
345 }
346 
347 /// Returns an array with all bins' minimum edges
348 /// The edges are arranges as xmin_1,ymin_1, xmin_2,ymin_2,....xmin_{nbin},ymin_{nbin}
350  if (fDataBins)
351  return &fBinMinEdges[0];
352  this->Warning("GetBinsMinEdges", "Binning kd-tree is nil. No bin edges retrieved.");
353  this->Info("GetBinsMinEdges", "Returning null pointer.");
354  return (Double_t*)0;
355 }
356 
357 /// Returns an array with all bins' maximum edges
358 /// The edges are arranges as xmax_1,ymax_1, xmax_2,ymax_2,....xmax_{nbin},ymax_{nbin}
360  // Returns the bins' maximum edges
361  if (fDataBins)
362  return &fBinMaxEdges[0];
363  this->Warning("GetBinsMaxEdges", "Binning kd-tree is nil. No bin edges retrieved.");
364  this->Info("GetBinsMaxEdges", "Returning null pointer.");
365  return (Double_t*)0;
366 }
367 
368 /// Returns a pair of an array with all bins minimum and maximum edges
369 std::pair<const Double_t*, const Double_t*> TKDTreeBinning::GetBinsEdges() const {
370  // Returns the bins' edges
371  if (fDataBins)
372  return std::make_pair(GetBinsMinEdges(), GetBinsMaxEdges());
373  this->Warning("GetBinsEdges", "Binning kd-tree is nil. No bin edges retrieved.");
374  this->Info("GetBinsEdges", "Returning null pointer pair.");
375  return std::make_pair((Double_t*)0, (Double_t*)0);
376 }
377 
378 /// Returns the bin's minimum edges. 'bin' is between 0 and fNBins - 1
380  if (fDataBins)
381  if (bin < fNBins)
382  return &fBinMinEdges[bin * fDim];
383  else
384  this->Warning("GetBinMinEdges", "No such bin. 'bin' is between 0 and %d", fNBins - 1);
385  else
386  this->Warning("GetBinMinEdges", "Binning kd-tree is nil. No bin edges retrieved.");
387  this->Info("GetBinMinEdges", "Returning null pointer.");
388  return (Double_t*)0;
389 }
390 
391 /// Returns the bin's maximum edges. 'bin' is between 0 and fNBins - 1
393  if (fDataBins)
394  if (bin < fNBins)
395  return &fBinMaxEdges[bin * fDim];
396  else
397  this->Warning("GetBinMaxEdges", "No such bin. 'bin' is between 0 and %d", fNBins - 1);
398  else
399  this->Warning("GetBinMaxEdges", "Binning kd-tree is nil. No bin edges retrieved.");
400  this->Info("GetBinMaxEdges", "Returning null pointer.");
401  return (Double_t*)0;
402 }
403 
404 /// Returns a pir with the bin's edges. 'bin' is between 0 and fNBins - 1
405 std::pair<const Double_t*, const Double_t*> TKDTreeBinning::GetBinEdges(UInt_t bin) const {
406  if (fDataBins)
407  if (bin < fNBins)
408  return std::make_pair(GetBinMinEdges(bin), GetBinMaxEdges(bin));
409  else
410  this->Warning("GetBinEdges", "No such bin. 'bin' is between 0 and %d", fNBins - 1);
411  else
412  this->Warning("GetBinEdges", "Binning kd-tree is nil. No bin edges retrieved.");
413  this->Info("GetBinEdges", "Returning null pointer pair.");
414  return std::make_pair((Double_t*)0, (Double_t*)0);
415 }
416 
417 /// Returns the number of bins
419  return fNBins;
420 }
421 
422 /// Returns the number of dimensions
424  return fDim;
425 }
426 
427 /// Returns the number of points in bin. 'bin' is between 0 and fNBins - 1
429  if(bin <= fNBins - 1)
430  return fBinsContent[bin];
431  this->Warning("GetBinContent", "No such bin. Returning 0.");
432  this->Info("GetBinContent", "'bin' is between 0 and %d.", fNBins - 1);
433  return 0;
434 }
435 
436 
437 /// Returns the kD-Tree structure of the binning
439  if (fDataBins)
440  return fDataBins;
441  this->Warning("GetTree", "Binning kd-tree is nil. No embedded kd-tree retrieved. Returning null pointer.");
442  return (TKDTreeID*)0;
443 }
444 
445 // Returns the data array in the dim coordinate. 'dim' is between 0 and fDim - 1
447  if(dim < fDim)
448  return &fData[dim*fDataSize];
449  this->Warning("GetDimData", "No such dimensional coordinate. No coordinate data retrieved. Returning null pointer.");
450  this->Info("GetDimData", "'dim' is between 0 and %d.", fDim - 1);
451  return 0;
452 }
453 
454 /// Returns the minimum of the data in the dim coordinate. 'dim' is between 0 and fDim - 1
456  if(dim < fDim)
457  return fDataThresholds[dim].first;
458  this->Warning("GetDataMin", "No such dimensional coordinate. No coordinate data minimum retrieved. Returning +inf.");
459  this->Info("GetDataMin", "'dim' is between 0 and %d.", fDim - 1);
460  return std::numeric_limits<Double_t>::infinity();
461 }
462 
463 /// Returns the maximum of the data in the dim coordinate. 'dim' is between 0 and fDim - 1
465  if(dim < fDim)
466  return fDataThresholds[dim].second;
467  this->Warning("GetDataMax", "No such dimensional coordinate. No coordinate data maximum retrieved. Returning -inf.");
468  this->Info("GetDataMax", "'dim' is between 0 and %d.", fDim - 1);
469  return -1 * std::numeric_limits<Double_t>::infinity();
470 }
471 
472 /// Returns the density in bin. 'bin' is between 0 and fNBins - 1
473 /// The density is the bin content/ bin volume
475  if(bin < fNBins) {
476  Double_t volume = GetBinVolume(bin);
477  if (!volume)
478  this->Warning("GetBinDensity", "Volume is null. Returning -1.");
479  return GetBinContent(bin) / volume;
480  }
481  this->Warning("GetBinDensity", "No such bin. Returning -1.");
482  this->Info("GetBinDensity", "'bin' is between 0 and %d.", fNBins - 1);
483  return -1.;
484 }
485 
486 /// Returns the (hyper)volume of bin. 'bin' is between 0 and fNBins - 1
488  if(bin < fNBins) {
489  std::pair<const Double_t*, const Double_t*> binEdges = GetBinEdges(bin);
490  Double_t volume = 1.;
491  for (UInt_t i = 0; i < fDim; ++i) {
492  volume *= (binEdges.second[i] - binEdges.first[i]);
493  }
494  return volume;
495  }
496  this->Warning("GetBinVolume", "No such bin. Returning 0.");
497  this->Info("GetBinVolume", "'bin' is between 0 and %d.", fNBins - 1);
498  return 0.;
499 }
500 
501 /// Returns a pointer to the vector of the bin edges for one dimensional binning only.
502 /// size of the vector is fNBins + 1 is the vector has been sorted in increasing bin edges
503 /// N.B : if one does not call SortOneDimBinEdges the bins are not ordered
504 const double * TKDTreeBinning::GetOneDimBinEdges() const {
505  if (fDim == 1) {
506  // no need to sort here because vector is already sorted in one dim
507  return &fBinMinEdges.front();
508  }
509  this->Warning("GetOneDimBinEdges", "Data is multidimensional. No sorted bin edges retrieved. Returning null pointer.");
510  this->Info("GetOneDimBinEdges", "This method can only be invoked if the data is a one dimensional set");
511  return 0;
512 }
513 
514 /// Sort the one-dimensional bin edges and retuns a pointer to them
516  if (fDim != 1) {
517  this->Warning("SortOneDimBinEdges", "Data is multidimensional. Cannot sorted bin edges. Returning null pointer.");
518  this->Info("SortOneDimBinEdges", "This method can only be invoked if the data is a one dimensional set");
519  return 0;
520  }
521  // order bins by increasing (or decreasing ) x positions
522  std::vector<UInt_t> indices(fNBins);
523  TMath::Sort( fNBins, &fBinMinEdges[0], &indices[0], !sortAsc );
524 
525  std::vector<Double_t> binMinEdges(fNBins );
526  std::vector<Double_t> binMaxEdges(fNBins );
527  std::vector<UInt_t> binContent(fNBins ); // reajust also content (not needed but better in general!)
528  fIndices.resize( fNBins );
529  for (UInt_t i = 0; i < fNBins; ++i) {
530  binMinEdges[i ] = fBinMinEdges[indices[i] ];
531  binMaxEdges[i ] = fBinMaxEdges[indices[i] ];
532  binContent[i] = fBinsContent[indices[i] ];
533  fIndices[indices[i] ] = i; // for the inverse transformation
534  }
535  fBinMinEdges.swap(binMinEdges);
536  fBinMaxEdges.swap(binMaxEdges);
537  fBinsContent.swap(binContent);
538 
539  fIsSorted = kTRUE;
540 
541  // Add also the upper(lower) edge to the min (max) list
542  if (sortAsc) {
543  fBinMinEdges.push_back(fBinMaxEdges.back());
544  fIsSortedAsc = kTRUE;
545  return &fBinMinEdges[0];
546  }
547  fBinMaxEdges.push_back(fBinMinEdges.back());
548  return &fBinMaxEdges[0];
549 
550 }
551 
552 /// Returns the geometric center of of the bin. 'bin' is between 0 and fNBins - 1
554  if(bin < fNBins) {
555  Double_t* result = new Double_t[fDim];
556  std::pair<const Double_t*, const Double_t*> binEdges = GetBinEdges(bin);
557  for (UInt_t i = 0; i < fDim; ++i) {
558  result[i] = (binEdges.second[i] + binEdges.first[i]) / 2.;
559  }
560  return result;
561  }
562  this->Warning("GetBinCenter", "No such bin. Returning null pointer.");
563  this->Info("GetBinCenter", "'bin' is between 0 and %d.", fNBins - 1);
564  return 0;
565 }
566 
567 /// Returns a pointer to the vector of the bin widths. 'bin' is between 0 and fNBins - 1
569  if(bin < fNBins) {
570  Double_t* result = new Double_t[fDim];
571  std::pair<const Double_t*, const Double_t*> binEdges = GetBinEdges(bin);
572  for (UInt_t i = 0; i < fDim; ++i) {
573  result[i] = (binEdges.second[i] - binEdges.first[i]);
574  }
575  return result;
576  }
577  this->Warning("GetBinWidth", "No such bin. Returning null pointer.");
578  this->Info("GetBinWidth", "'bin' is between 0 and %d.", fNBins - 1);
579  return 0;
580 }
581 
582 /// Return the bin with maximum density
584  if (fIsSorted) {
585  if (fIsSortedAsc)
586  return fNBins - 1;
587  else return 0;
588  }
589  UInt_t* indices = new UInt_t[fNBins];
590  for (UInt_t i = 0; i < fNBins; ++i)
591  indices[i] = i;
592  UInt_t result = *std::max_element(indices, indices + fNBins, CompareAsc(this));
593  delete [] indices;
594  return result;
595 }
596 
597 /// Return the bin with minimum density
599  if (fIsSorted) {
600  if (!fIsSortedAsc)
601  return fNBins - 1;
602  else return 0;
603  }
604  UInt_t* indices = new UInt_t[fNBins];
605  for (UInt_t i = 0; i < fNBins; ++i)
606  indices[i] = i;
607  UInt_t result = *std::min_element(indices, indices + fNBins, CompareAsc(this));
608  delete [] indices;
609  return result;
610 }
611 
612 /// Fill the bin data set (class ROOT::Fit::BinData) with the result of the TKDTree binning
614  if (!fDataBins) return;
615  data.Initialize(fNBins, fDim);
616  for (unsigned int i = 0; i < fNBins; ++i) {
617  data.Add( GetBinMinEdges(i), GetBinDensity(i), std::sqrt(double(GetBinContent(i) ))/ GetBinVolume(i) );
618  data.AddBinUpEdge(GetBinMaxEdges(i) );
619  }
620 }
621 
622 /// find the corresponding bin index given the coordinate of a point
624 
625  Int_t inode = fDataBins->FindNode(point);
626  // find node return the index in the total nodes and the bins are only the terminal ones
627  // so we subtract all the non-terminal nodes
628  inode -= fDataBins->GetNNodes();
629  R__ASSERT( inode >= 0);
630  UInt_t bin = inode;
631 
632  if (!fIsSorted) return bin;
633  //return std::distance(fIndices.begin(), std::find(fIndices.begin(), fIndices.end(), bin ) );
634  return fIndices[bin];
635 }
636 
637 /// Return the corresponding point belonging to the bin i
638 std::vector<std::vector<Double_t> > TKDTreeBinning::GetPointsInBin(UInt_t bin) const {
639  std::vector<Double_t> point(fDim);
640  std::vector< std::vector<Double_t> > thePoints;
641  if (fData.size() == 0) {
642  Error("GetPointsInBin","Internal data set is not valid");
643  return thePoints;
644  }
645  if (!fDataBins) {
646  Error("GetPointsInBin","Internal TKDTree is not valid");
647  return thePoints;
648  }
649  if (bin >= fNBins) {
650  Error("GetPointsInBin","Invalid bin number");
651  return thePoints;
652  }
653  // need to find the bin number including the non-terminal node
654  int inode = bin + fDataBins->GetNNodes();
655  auto indices = fDataBins->GetPointsIndexes(inode);
656  int npoints = fDataBins->GetNPointsNode(inode);
657  thePoints.resize(npoints);
658  for (int ipoint = 0; ipoint < npoints; ++ipoint) {
659  for (unsigned int idim = 0; idim < fDim; ++idim) {
660  point[idim] = fData[idim*fDataSize+indices[ipoint] ];
661  }
662  thePoints[ipoint] = point;
663  }
664  return thePoints;
665 }
666 
667 ////////////////////////////////////////////////////////////////////////////////
668 /// Stream a class object.
669 void TKDTreeBinning::Streamer(TBuffer &b) {
670  if (b.IsReading() ) {
671  UInt_t R__s, R__c;
672  Version_t v = b.ReadVersion(&R__s, &R__c);
673  b.ReadClassBuffer(TKDTreeBinning::Class(), this, v, R__s, R__c);
674  // we need to rebuild the TKDTree
675  if (fDataBins) delete fDataBins;
676  SetNBins(fNBins);
677  }
678  else {
679  // case of writing
681  //std::cout << "writing npar = " << GetNpar() << std::endl;
682  }
683 }
TKDTreeID * GetTree() const
Returns the kD-Tree structure of the binning.
std::vector< std::vector< std::pair< Bool_t, Bool_t > > > fCheckedBinEdges
Minimum and maximum data values.
Bool_t IsReading() const
Definition: TBuffer.h:83
void Initialize(unsigned int maxpoints, unsigned int dim=1, ErrorType err=kValueError)
preallocate a data set with given size , dimension and error type (to get the full point size) If the...
Definition: BinData.cxx:216
std::vector< Double_t > fBinMaxEdges
The minimum values for the bins&#39; edges for each dimension.
virtual Int_t WriteClassBuffer(const TClass *cl, void *pointer)=0
virtual void Info(const char *method, const char *msgfmt,...) const
Issue info message.
Definition: TObject.cxx:899
void SetData(Index npoints, Index ndim, UInt_t bsize, Value **data)
Set the data array. See the constructor function comments for details.
Definition: TKDTree.cxx:923
const Double_t * GetBinsMinEdges() const
Returns an array with all bins&#39; minimum edges The edges are arranges as xmin_1,ymin_1, xmin_2,ymin_2,....xmin_{nbin},ymin_{nbin}.
Index FindNode(const Value *point) const
returns the index of the terminal node to which point belongs (index in the fAxis, fValue, etc arrays) returns -1 in case of failure
Definition: TKDTree.cxx:677
void SetBinMinMaxEdges(Double_t *binEdges)
Index * GetPointsIndexes(Int_t node) const
return the indices of the points in that terminal node for all the nodes except last, the size is fBucketSize for the last node it&#39;s fOffsetfBucketSize
Definition: TKDTree.cxx:817
Double_t GetDataMax(UInt_t dim) const
Returns the maximum of the data in the dim coordinate. &#39;dim&#39; is between 0 and fDim - 1...
short Version_t
Definition: RtypesCore.h:61
friend struct CompareDesc
! Predicate for descending sort
std::vector< UInt_t > fBinsContent
Flags if the bin edges are sorted densitywise (or by bin-edge for 1D) in ascending order...
Bool_t TestBit(UInt_t f) const
Definition: TObject.h:157
UInt_t GetBinMinDensity() const
Return the bin with minimum density.
void SetData(Double_t *data)
Disallowed assign operator.
const Double_t * GetBinsMaxEdges() const
Returns an array with all bins&#39; maximum edges The edges are arranges as xmax_1,ymax_1, xmax_2,ymax_2,....xmax_{nbin},ymax_{nbin}.
std::pair< const Double_t *, const Double_t * > GetBinEdges(UInt_t bin) const
Returns a pir with the bin&#39;s edges. &#39;bin&#39; is between 0 and fNBins - 1.
Buffer base class used for serializing objects.
Definition: TBuffer.h:42
#define R__ASSERT(e)
Definition: TError.h:98
UInt_t GetNBins() const
Returns the number of bins.
TKDTreeBinning()
Default constructor (for I/O)
int Int_t
Definition: RtypesCore.h:41
bool Bool_t
Definition: RtypesCore.h:59
const Bool_t kFALSE
Definition: Rtypes.h:92
void SetNBins(UInt_t bins)
Sets binning inner structure.
Index GetNPointsNode(Int_t node) const
Get number of points in this node for all the terminal nodes except last, the size is fBucketSize for...
Definition: TKDTree.cxx:900
STL namespace.
const char * Class
Definition: TXMLSetup.cxx:64
void AddBinUpEdge(const double *xup)
add the bin width data, a pointer to an array with the bin upper edge information.
Definition: BinData.cxx:452
void Build()
Build the kd-tree.
Definition: TKDTree.cxx:411
void SetBit(UInt_t f, Bool_t set)
Set or unset the user status bits as specified in f.
Definition: TObject.cxx:739
UInt_t fDataSize
The data dimension.
Bool_t fIsSortedAsc
Flags if the bin edges are sorted densitywise (or by bin endges in case of 1-dim ) ...
double sqrt(double)
Index GetBucketSize()
Definition: TKDTree.h:50
~TKDTreeBinning()
Class&#39;s destructor.
void Sort(Index n, const Element *a, Index *index, Bool_t down=kTRUE)
Definition: TMath.h:989
const Double_t * SortOneDimBinEdges(Bool_t sortAsc=kTRUE)
Sort the one-dimensional bin edges and retuns a pointer to them.
<- TKDTreeBinning - A class providing multidimensional binning ->
Double_t GetDataMin(UInt_t dim) const
Returns the minimum of the data in the dim coordinate. &#39;dim&#39; is between 0 and fDim - 1...
Class implementing a kd-tree.
Definition: TKDTree.h:11
TKDTreeID * fDataBins
Index of the bins in the kd-tree (needed when bins are sorted)
const Double_t * GetBinCenter(UInt_t bin) const
Returns the geometric center of of the bin. &#39;bin&#39; is between 0 and fNBins - 1.
void ReadjustMaxBinEdges(Double_t *binEdges)
std::vector< Double_t > fBinMinEdges
[fDataSize*fDim] The data from which a KDTree partition is computed for binning
Double_t GetBinDensity(UInt_t bin) const
Returns the density in bin.
Int_t GetNNodes() const
Definition: TKDTree.h:35
VecExpr< UnaryOp< Fabs< T >, VecExpr< A, T, D >, T >, T, D > fabs(const VecExpr< A, T, D > &rhs)
std::vector< Double_t > fData
UInt_t GetBinContent(UInt_t bin) const
Returns the number of points in bin. &#39;bin&#39; is between 0 and fNBins - 1.
SVector< double, 2 > v
Definition: Dict.h:5
std::vector< UInt_t > fIndices
The maximum values for the bins&#39; edges for each dimension.
const Double_t * GetBinMinEdges(UInt_t bin) const
Returns the bin&#39;s minimum edges. &#39;bin&#39; is between 0 and fNBins - 1.
unsigned int UInt_t
Definition: RtypesCore.h:42
virtual void Error(const char *method, const char *msgfmt,...) const
Issue error message.
Definition: TObject.cxx:925
const Double_t * GetBinWidth(UInt_t bin) const
Returns a pointer to the vector of the bin widths. &#39;bin&#39; is between 0 and fNBins - 1...
friend struct CompareAsc
! Predicate for ascending sort
const Double_t * GetDimData(UInt_t dim) const
Class describing the binned data sets : vectors of x coordinates, y values and optionally error on y ...
Definition: BinData.h:61
REAL epsilon
Definition: triangle.c:617
std::vector< std::pair< Double_t, Double_t > > fDataThresholds
The data size.
TKDTree< Int_t, Double_t > TKDTreeID
Definition: TKDTree.h:106
void SortBinsByDensity(Bool_t sortAsc=kTRUE)
Sorts bins by their density.
virtual Int_t ReadClassBuffer(const TClass *cl, void *pointer, const TClass *onfile_class=0)=0
#define ClassImp(name)
Definition: Rtypes.h:279
void FillBinData(ROOT::Fit::BinData &data) const
Fill the bin data set (class ROOT::Fit::BinData) with the result of the TKDTree binning.
double Double_t
Definition: RtypesCore.h:55
void Add(double x, double y)
add one dim data with only coordinate and values
Definition: BinData.cxx:265
UInt_t GetBinMaxDensity() const
Return the bin with maximum density.
std::pair< const Double_t *, const Double_t * > GetBinsEdges() const
Returns a pair of an array with all bins minimum and maximum edges.
void SetCommonBinEdges(Double_t *binEdges)
const Double_t * GetBinMaxEdges(UInt_t bin) const
Returns the bin&#39;s maximum edges. &#39;bin&#39; is between 0 and fNBins - 1.
void ReadjustMinBinEdges(Double_t *binEdges)
UInt_t fDim
The number of bins.
#define nullptr
Definition: Rtypes.h:87
std::vector< std::map< Double_t, std::vector< UInt_t > > > fCommonBinEdges
! Auxiliary structure for readjusting the bin edges. Keeps the common bin boundaries ...
you should not use this method at all Int_t Int_t Double_t Double_t Double_t Int_t Double_t Double_t Double_t Double_t b
Definition: TRolke.cxx:630
Double_t GetBinVolume(UInt_t bin) const
Returns the (hyper)volume of bin. &#39;bin&#39; is between 0 and fNBins - 1.
double result[121]
Definition: first.py:1
UInt_t FindBin(const Double_t *point) const
find the corresponding bin index given the coordinate of a point
UInt_t GetDim() const
Returns the number of dimensions.
const Bool_t kTRUE
Definition: Rtypes.h:91
std::vector< std::vector< Double_t > > GetPointsInBin(UInt_t bin) const
Return the corresponding point belonging to the bin i.
Value * GetBoundary(const Int_t node)
Get a boundary.
Definition: TKDTree.cxx:1211
virtual void Warning(const char *method, const char *msgfmt,...) const
Issue warning message.
Definition: TObject.cxx:911
virtual Version_t ReadVersion(UInt_t *start=0, UInt_t *bcnt=0, const TClass *cl=0)=0
const Double_t * GetOneDimBinEdges() const
Returns a pointer to the vector of the bin edges for one dimensional binning only.