Logo ROOT   6.10/09
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), fNBins (nBins), 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),
115  fDataBins(nullptr),
116  fNBins (0),
117  fDim(0),
118  fDataSize(0),
120 {}
121 
122 /// Class's destructor
124  // if (fData) delete[] fData;
125  if (fDataBins) delete fDataBins;
126 }
127 
128 /// Sets binning inner structure
130  fNBins = bins;
131  if (fDim && fNBins && fDataSize) {
132  if (fDataSize / fNBins) {
133  Bool_t remainingData = fDataSize % fNBins;
134  if (remainingData) {
135  fNBins += 1;
136  this->Info("SetNBins", "Number of bins is not enough to hold the data. Extra bin added.");
137  }
138  fDataBins = new TKDTreeID(fDataSize, fDim, fDataSize / (fNBins - remainingData)); // TKDTree input is data size, data dimension and the content size of bins ("bucket" size)
139  SetTreeData();
140  fDataBins->Build();
141  SetBinsEdges();
142  SetBinsContent();
143  } else {
144  fDataBins = (TKDTreeID*)0;
145  this->Warning("SetNBins", "Number of bins is bigger than data size. Nothing is built.");
146  }
147  } else {
148  fDataBins = (TKDTreeID*)0;
149  if (!fDim)
150  this->Warning("SetNBins", "Data dimension is nil. Nothing is built.");
151  if (!fNBins)
152  this->Warning("SetNBins", "Number of bins is nil. Nothing is built.");
153  if (!fDataSize)
154  this->Warning("SetNBins", "Data size is nil. Nothing is built.");
155  }
156 }
157 
158 /// Sorts bins by their density
160  if (fDim == 1) {
161  // in one dim they are already sorted (no need to do anything)
162  return;
163  } else {
164  std::vector<UInt_t> indices(fNBins); // vector for indices (for the inverse transformation)
165  for (UInt_t i = 0; i < fNBins; ++i)
166  indices[i] = i;
167  if (sortAsc) {
168  std::sort(indices.begin(), indices.end(), CompareAsc(this));
170  } else {
171  std::sort(indices.begin(), indices.end(), CompareDesc(this));
173  }
174  std::vector<Double_t> binMinEdges(fNBins * fDim);
175  std::vector<Double_t> binMaxEdges(fNBins * fDim);
176  std::vector<UInt_t> binContent(fNBins ); // reajust also content (not needed bbut better in general!)
177  fIndices.resize(fNBins);
178  for (UInt_t i = 0; i < fNBins; ++i) {
179  for (UInt_t j = 0; j < fDim; ++j) {
180  binMinEdges[i * fDim + j] = fBinMinEdges[indices[i] * fDim + j];
181  binMaxEdges[i * fDim + j] = fBinMaxEdges[indices[i] * fDim + j];
182  }
183  binContent[i] = fBinsContent[indices[i]];
184  fIndices[indices[i]] = i;
185  }
186  fBinMinEdges.swap(binMinEdges);
187  fBinMaxEdges.swap(binMaxEdges);
188  fBinsContent.swap(binContent);
189 
190  // not needed anymore if readjusting bin content all the time
191  // re-adjust content of extra bins if exists
192  // since it is different than the others
193  // if ( fDataSize % fNBins != 0) {
194  // UInt_t k = 0;
195  // Bool_t found = kFALSE;
196  // while (!found) {
197  // if (indices[k] == fNBins - 1) {
198  // found = kTRUE;
199  // break;
200  // }
201  // ++k;
202  // }
203  // fBinsContent[fNBins - 1] = fDataBins->GetBucketSize();
204  // fBinsContent[k] = fDataSize % fNBins-1;
205  // }
206 
207  fIsSorted = kTRUE;
208  }
209 }
210 
212  // Sets the data and finds minimum and maximum by dimensional coordinate
213  fData.resize(fDim*fDataSize);
214  auto first = fData.begin();
215  for (UInt_t i = 0; i < fDim; ++i) {
216  for (UInt_t j = 0; j < fDataSize; ++j) {
217  fData[i*fDataSize+j] = data[i * fDataSize + j];
218  }
219  auto end = first+fDataSize;
220  fDataThresholds[i] = std::make_pair(*std::min_element(first, end), *std::max_element(first,end));
221  first = end;
222  }
223 }
224 void TKDTreeBinning::SetData(const std::vector<double>& data) {
225  // Sets the data and finds minimum and maximum by dimensional coordinate
226  fData = data;
227  auto first = fData.begin();
228  // find min/max
229  for (UInt_t i = 0; i < fDim; ++i) {
230  auto end = first+fDataSize;
231  fDataThresholds[i] = std::make_pair(*std::min_element(first, end), *std::max_element(first,end));
232  first = end;
233  }
234 }
235 
237  // Sets the data for constructing the kD-tree
238  for (UInt_t i = 0; i < fDim; ++i)
239  fDataBins->SetData(i, &fData[i*fDataSize]);
240 }
241 
243  // Sets the bins' content
244  fBinsContent.reserve(fNBins);
245  for (UInt_t i = 0; i < fNBins; ++i)
247  if ( fDataSize % fNBins != 0 )
248  fBinsContent[fNBins - 1] = fDataSize % (fNBins-1);
249 }
250 
252  // Sets the bins' edges
253  //Double_t* rawBinEdges = fDataBins->GetBoundaryExact(fDataBins->GetNNodes());
254  Double_t* rawBinEdges = fDataBins->GetBoundary(fDataBins->GetNNodes());
255  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)));
256  fCommonBinEdges = std::vector<std::map<Double_t, std::vector<UInt_t> > >(fDim, std::map<Double_t, std::vector<UInt_t> >());
257  SetCommonBinEdges(rawBinEdges);
258  if (TestBit(kAdjustBinEdges) ) {
259  ReadjustMinBinEdges(rawBinEdges);
260  ReadjustMaxBinEdges(rawBinEdges);
261  }
262  SetBinMinMaxEdges(rawBinEdges);
263  fCommonBinEdges.clear();
264  fCheckedBinEdges.clear();
265 }
266 
268  // Sets the bins' minimum and maximum edges
269  fBinMinEdges.reserve(fNBins * fDim);
270  fBinMaxEdges.reserve(fNBins * fDim);
271  for (UInt_t i = 0; i < fNBins; ++i) {
272  for (UInt_t j = 0; j < fDim; ++j) {
273  fBinMinEdges.push_back(binEdges[(i * fDim + j) * 2]);
274  fBinMaxEdges.push_back(binEdges[(i * fDim + j) * 2 + 1]);
275  }
276  }
277 }
278 
280  // Sets indexing on the bin edges which have common boundaries
281  for (UInt_t i = 0; i < fDim; ++i) {
282  for (UInt_t j = 0; j < fNBins; ++j) {
283  Double_t binEdge = binEdges[(j * fDim + i) * 2];
284  if(fCommonBinEdges[i].find(binEdge) == fCommonBinEdges[i].end()) {
285  std::vector<UInt_t> commonBinEdges;
286  for (UInt_t k = 0; k < fNBins; ++k) {
287  UInt_t minBinEdgePos = (k * fDim + i) * 2;
288  if (std::fabs(binEdge - binEdges[minBinEdgePos]) < std::numeric_limits<Double_t>::epsilon())
289  commonBinEdges.push_back(minBinEdgePos);
290  UInt_t maxBinEdgePos = ++minBinEdgePos;
291  if (std::fabs(binEdge - binEdges[maxBinEdgePos]) < std::numeric_limits<Double_t>::epsilon())
292  commonBinEdges.push_back(maxBinEdgePos);
293  }
294  fCommonBinEdges[i][binEdge] = commonBinEdges;
295  }
296  }
297  }
298 }
299 
301  // Readjusts the bins' minimum edge by shifting it slightly lower
302  // to avoid overlapping with the data
303  for (UInt_t i = 0; i < fDim; ++i) {
304  for (UInt_t j = 0; j < fNBins; ++j) {
305  if (!fCheckedBinEdges[i][j].first) {
306  Double_t binEdge = binEdges[(j * fDim + i) * 2];
307  Double_t adjustedBinEdge = binEdge;
308  double eps = -10*std::numeric_limits<Double_t>::epsilon();
309  if (adjustedBinEdge != 0)
310  adjustedBinEdge *= (1. + eps);
311  else
312  adjustedBinEdge += eps;
313 
314  for (UInt_t k = 0; k < fCommonBinEdges[i][binEdge].size(); ++k) {
315  UInt_t binEdgePos = fCommonBinEdges[i][binEdge][k];
316  Bool_t isMinBinEdge = binEdgePos % 2 == 0;
317  UInt_t bin = isMinBinEdge ? (binEdgePos / 2 - i) / fDim : ((binEdgePos - 1) / 2 - i) / fDim;
318  binEdges[binEdgePos] = adjustedBinEdge;
319  if (isMinBinEdge)
320  fCheckedBinEdges[i][bin].first = kTRUE;
321  else
322  fCheckedBinEdges[i][bin].second = kTRUE;
323  }
324  }
325  }
326  }
327 }
328 
330  // Readjusts the bins' maximum edge
331  // and shift it sligtly higher
332  for (UInt_t i = 0; i < fDim; ++i) {
333  for (UInt_t j = 0; j < fNBins; ++j) {
334  if (!fCheckedBinEdges[i][j].second) {
335  Double_t& binEdge = binEdges[(j * fDim + i) * 2 + 1];
336  double eps = 10*std::numeric_limits<Double_t>::epsilon();
337  if (binEdge != 0)
338  binEdge *= (1. + eps);
339  else
340  binEdge += eps;
341 
342 
343  }
344  }
345  }
346 }
347 
348 /// Returns an array with all bins' minimum edges
349 /// The edges are arranges as xmin_1,ymin_1, xmin_2,ymin_2,....xmin_{nbin},ymin_{nbin}
351  if (fDataBins)
352  return &fBinMinEdges[0];
353  this->Warning("GetBinsMinEdges", "Binning kd-tree is nil. No bin edges retrieved.");
354  this->Info("GetBinsMinEdges", "Returning null pointer.");
355  return (Double_t*)0;
356 }
357 
358 /// Returns an array with all bins' maximum edges
359 /// The edges are arranges as xmax_1,ymax_1, xmax_2,ymax_2,....xmax_{nbin},ymax_{nbin}
361  // Returns the bins' maximum edges
362  if (fDataBins)
363  return &fBinMaxEdges[0];
364  this->Warning("GetBinsMaxEdges", "Binning kd-tree is nil. No bin edges retrieved.");
365  this->Info("GetBinsMaxEdges", "Returning null pointer.");
366  return (Double_t*)0;
367 }
368 
369 /// Returns a pair of an array with all bins minimum and maximum edges
370 std::pair<const Double_t*, const Double_t*> TKDTreeBinning::GetBinsEdges() const {
371  // Returns the bins' edges
372  if (fDataBins)
373  return std::make_pair(GetBinsMinEdges(), GetBinsMaxEdges());
374  this->Warning("GetBinsEdges", "Binning kd-tree is nil. No bin edges retrieved.");
375  this->Info("GetBinsEdges", "Returning null pointer pair.");
376  return std::make_pair((Double_t*)0, (Double_t*)0);
377 }
378 
379 /// Returns the bin's minimum edges. 'bin' is between 0 and fNBins - 1
381  if (fDataBins)
382  if (bin < fNBins)
383  return &fBinMinEdges[bin * fDim];
384  else
385  this->Warning("GetBinMinEdges", "No such bin. 'bin' is between 0 and %d", fNBins - 1);
386  else
387  this->Warning("GetBinMinEdges", "Binning kd-tree is nil. No bin edges retrieved.");
388  this->Info("GetBinMinEdges", "Returning null pointer.");
389  return (Double_t*)0;
390 }
391 
392 /// Returns the bin's maximum edges. 'bin' is between 0 and fNBins - 1
394  if (fDataBins)
395  if (bin < fNBins)
396  return &fBinMaxEdges[bin * fDim];
397  else
398  this->Warning("GetBinMaxEdges", "No such bin. 'bin' is between 0 and %d", fNBins - 1);
399  else
400  this->Warning("GetBinMaxEdges", "Binning kd-tree is nil. No bin edges retrieved.");
401  this->Info("GetBinMaxEdges", "Returning null pointer.");
402  return (Double_t*)0;
403 }
404 
405 /// Returns a pir with the bin's edges. 'bin' is between 0 and fNBins - 1
406 std::pair<const Double_t*, const Double_t*> TKDTreeBinning::GetBinEdges(UInt_t bin) const {
407  if (fDataBins)
408  if (bin < fNBins)
409  return std::make_pair(GetBinMinEdges(bin), GetBinMaxEdges(bin));
410  else
411  this->Warning("GetBinEdges", "No such bin. 'bin' is between 0 and %d", fNBins - 1);
412  else
413  this->Warning("GetBinEdges", "Binning kd-tree is nil. No bin edges retrieved.");
414  this->Info("GetBinEdges", "Returning null pointer pair.");
415  return std::make_pair((Double_t*)0, (Double_t*)0);
416 }
417 
418 /// Returns the number of bins
420  return fNBins;
421 }
422 
423 /// Returns the number of dimensions
425  return fDim;
426 }
427 
428 /// Returns the number of points in bin. 'bin' is between 0 and fNBins - 1
430  if(bin <= fNBins - 1)
431  return fBinsContent[bin];
432  this->Warning("GetBinContent", "No such bin. Returning 0.");
433  this->Info("GetBinContent", "'bin' is between 0 and %d.", fNBins - 1);
434  return 0;
435 }
436 
437 
438 /// Returns the kD-Tree structure of the binning
440  if (fDataBins)
441  return fDataBins;
442  this->Warning("GetTree", "Binning kd-tree is nil. No embedded kd-tree retrieved. Returning null pointer.");
443  return (TKDTreeID*)0;
444 }
445 
446 // Returns the data array in the dim coordinate. 'dim' is between 0 and fDim - 1
448  if(dim < fDim)
449  return &fData[dim*fDataSize];
450  this->Warning("GetDimData", "No such dimensional coordinate. No coordinate data retrieved. Returning null pointer.");
451  this->Info("GetDimData", "'dim' is between 0 and %d.", fDim - 1);
452  return 0;
453 }
454 
455 /// Returns the minimum of the data in the dim coordinate. 'dim' is between 0 and fDim - 1
457  if(dim < fDim)
458  return fDataThresholds[dim].first;
459  this->Warning("GetDataMin", "No such dimensional coordinate. No coordinate data minimum retrieved. Returning +inf.");
460  this->Info("GetDataMin", "'dim' is between 0 and %d.", fDim - 1);
461  return std::numeric_limits<Double_t>::infinity();
462 }
463 
464 /// Returns the maximum of the data in the dim coordinate. 'dim' is between 0 and fDim - 1
466  if(dim < fDim)
467  return fDataThresholds[dim].second;
468  this->Warning("GetDataMax", "No such dimensional coordinate. No coordinate data maximum retrieved. Returning -inf.");
469  this->Info("GetDataMax", "'dim' is between 0 and %d.", fDim - 1);
470  return -1 * std::numeric_limits<Double_t>::infinity();
471 }
472 
473 /// Returns the density in bin. 'bin' is between 0 and fNBins - 1
474 /// The density is the bin content/ bin volume
476  if(bin < fNBins) {
477  Double_t volume = GetBinVolume(bin);
478  if (!volume)
479  this->Warning("GetBinDensity", "Volume is null. Returning -1.");
480  return GetBinContent(bin) / volume;
481  }
482  this->Warning("GetBinDensity", "No such bin. Returning -1.");
483  this->Info("GetBinDensity", "'bin' is between 0 and %d.", fNBins - 1);
484  return -1.;
485 }
486 
487 /// Returns the (hyper)volume of bin. 'bin' is between 0 and fNBins - 1
489  if(bin < fNBins) {
490  std::pair<const Double_t*, const Double_t*> binEdges = GetBinEdges(bin);
491  Double_t volume = 1.;
492  for (UInt_t i = 0; i < fDim; ++i) {
493  volume *= (binEdges.second[i] - binEdges.first[i]);
494  }
495  return volume;
496  }
497  this->Warning("GetBinVolume", "No such bin. Returning 0.");
498  this->Info("GetBinVolume", "'bin' is between 0 and %d.", fNBins - 1);
499  return 0.;
500 }
501 
502 /// Returns a pointer to the vector of the bin edges for one dimensional binning only.
503 /// size of the vector is fNBins + 1 is the vector has been sorted in increasing bin edges
504 /// N.B : if one does not call SortOneDimBinEdges the bins are not ordered
505 const double * TKDTreeBinning::GetOneDimBinEdges() const {
506  if (fDim == 1) {
507  // no need to sort here because vector is already sorted in one dim
508  return &fBinMinEdges.front();
509  }
510  this->Warning("GetOneDimBinEdges", "Data is multidimensional. No sorted bin edges retrieved. Returning null pointer.");
511  this->Info("GetOneDimBinEdges", "This method can only be invoked if the data is a one dimensional set");
512  return 0;
513 }
514 
515 /// Sort the one-dimensional bin edges and retuns a pointer to them
517  if (fDim != 1) {
518  this->Warning("SortOneDimBinEdges", "Data is multidimensional. Cannot sorted bin edges. Returning null pointer.");
519  this->Info("SortOneDimBinEdges", "This method can only be invoked if the data is a one dimensional set");
520  return 0;
521  }
522  // order bins by increasing (or decreasing ) x positions
523  std::vector<UInt_t> indices(fNBins);
524  TMath::Sort( fNBins, &fBinMinEdges[0], &indices[0], !sortAsc );
525 
526  std::vector<Double_t> binMinEdges(fNBins );
527  std::vector<Double_t> binMaxEdges(fNBins );
528  std::vector<UInt_t> binContent(fNBins ); // reajust also content (not needed but better in general!)
529  fIndices.resize( fNBins );
530  for (UInt_t i = 0; i < fNBins; ++i) {
531  binMinEdges[i ] = fBinMinEdges[indices[i] ];
532  binMaxEdges[i ] = fBinMaxEdges[indices[i] ];
533  binContent[i] = fBinsContent[indices[i] ];
534  fIndices[indices[i] ] = i; // for the inverse transformation
535  }
536  fBinMinEdges.swap(binMinEdges);
537  fBinMaxEdges.swap(binMaxEdges);
538  fBinsContent.swap(binContent);
539 
540  fIsSorted = kTRUE;
541 
542  // Add also the upper(lower) edge to the min (max) list
543  if (sortAsc) {
544  fBinMinEdges.push_back(fBinMaxEdges.back());
545  fIsSortedAsc = kTRUE;
546  return &fBinMinEdges[0];
547  }
548  fBinMaxEdges.push_back(fBinMinEdges.back());
549  return &fBinMaxEdges[0];
550 
551 }
552 
553 /// Returns the geometric center of of the bin. 'bin' is between 0 and fNBins - 1
555  if(bin < fNBins) {
556  Double_t* result = new Double_t[fDim];
557  std::pair<const Double_t*, const Double_t*> binEdges = GetBinEdges(bin);
558  for (UInt_t i = 0; i < fDim; ++i) {
559  result[i] = (binEdges.second[i] + binEdges.first[i]) / 2.;
560  }
561  return result;
562  }
563  this->Warning("GetBinCenter", "No such bin. Returning null pointer.");
564  this->Info("GetBinCenter", "'bin' is between 0 and %d.", fNBins - 1);
565  return 0;
566 }
567 
568 /// Returns a pointer to the vector of the bin widths. 'bin' is between 0 and fNBins - 1
570  if(bin < fNBins) {
571  Double_t* result = new Double_t[fDim];
572  std::pair<const Double_t*, const Double_t*> binEdges = GetBinEdges(bin);
573  for (UInt_t i = 0; i < fDim; ++i) {
574  result[i] = (binEdges.second[i] - binEdges.first[i]);
575  }
576  return result;
577  }
578  this->Warning("GetBinWidth", "No such bin. Returning null pointer.");
579  this->Info("GetBinWidth", "'bin' is between 0 and %d.", fNBins - 1);
580  return 0;
581 }
582 
583 /// Return the bin with maximum density
585  if (fIsSorted) {
586  if (fIsSortedAsc)
587  return fNBins - 1;
588  else return 0;
589  }
590  UInt_t* indices = new UInt_t[fNBins];
591  for (UInt_t i = 0; i < fNBins; ++i)
592  indices[i] = i;
593  UInt_t result = *std::max_element(indices, indices + fNBins, CompareAsc(this));
594  delete [] indices;
595  return result;
596 }
597 
598 /// Return the bin with minimum density
600  if (fIsSorted) {
601  if (!fIsSortedAsc)
602  return fNBins - 1;
603  else return 0;
604  }
605  UInt_t* indices = new UInt_t[fNBins];
606  for (UInt_t i = 0; i < fNBins; ++i)
607  indices[i] = i;
608  UInt_t result = *std::min_element(indices, indices + fNBins, CompareAsc(this));
609  delete [] indices;
610  return result;
611 }
612 
613 /// Fill the bin data set (class ROOT::Fit::BinData) with the result of the TKDTree binning
615  if (!fDataBins) return;
616  data.Initialize(fNBins, fDim);
617  for (unsigned int i = 0; i < fNBins; ++i) {
618  data.Add( GetBinMinEdges(i), GetBinDensity(i), std::sqrt(double(GetBinContent(i) ))/ GetBinVolume(i) );
619  data.AddBinUpEdge(GetBinMaxEdges(i) );
620  }
621 }
622 
623 /// find the corresponding bin index given the coordinate of a point
625 
626  Int_t inode = fDataBins->FindNode(point);
627  // find node return the index in the total nodes and the bins are only the terminal ones
628  // so we subtract all the non-terminal nodes
629  inode -= fDataBins->GetNNodes();
630  R__ASSERT( inode >= 0);
631  UInt_t bin = inode;
632 
633  if (!fIsSorted) return bin;
634  //return std::distance(fIndices.begin(), std::find(fIndices.begin(), fIndices.end(), bin ) );
635  return fIndices[bin];
636 }
637 
638 /// Return the corresponding point belonging to the bin i
639 std::vector<std::vector<Double_t> > TKDTreeBinning::GetPointsInBin(UInt_t bin) const {
640  std::vector<Double_t> point(fDim);
641  std::vector< std::vector<Double_t> > thePoints;
642  if (fData.size() == 0) {
643  Error("GetPointsInBin","Internal data set is not valid");
644  return thePoints;
645  }
646  if (!fDataBins) {
647  Error("GetPointsInBin","Internal TKDTree is not valid");
648  return thePoints;
649  }
650  if (bin >= fNBins) {
651  Error("GetPointsInBin","Invalid bin number");
652  return thePoints;
653  }
654  // need to find the bin number including the non-terminal node
655  int inode = bin + fDataBins->GetNNodes();
656  auto indices = fDataBins->GetPointsIndexes(inode);
657  int npoints = fDataBins->GetNPointsNode(inode);
658  thePoints.resize(npoints);
659  for (int ipoint = 0; ipoint < npoints; ++ipoint) {
660  for (unsigned int idim = 0; idim < fDim; ++idim) {
661  point[idim] = fData[idim*fDataSize+indices[ipoint] ];
662  }
663  thePoints[ipoint] = point;
664  }
665  return thePoints;
666 }
667 
668 ////////////////////////////////////////////////////////////////////////////////
669 /// Stream a class object.
670 void TKDTreeBinning::Streamer(TBuffer &b) {
671  if (b.IsReading() ) {
672  UInt_t R__s, R__c;
673  Version_t v = b.ReadVersion(&R__s, &R__c);
674  b.ReadClassBuffer(TKDTreeBinning::Class(), this, v, R__s, R__c);
675  // we need to rebuild the TKDTree
676  if (fDataBins) delete fDataBins;
677  SetNBins(fNBins);
678  }
679  else {
680  // case of writing
682  //std::cout << "writing npar = " << GetNpar() << std::endl;
683  }
684 }
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:81
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:847
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:159
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:40
#define R__ASSERT(e)
Definition: TError.h:96
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
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.
void AddBinUpEdge(const double *xup)
add the bin width data, a pointer to an array with the bin upper edge information.
Definition: BinData.cxx:596
TRObject operator()(const T1 &t1) const
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:687
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)
void Class()
Definition: Class.C:29
Index GetBucketSize()
Definition: TKDTree.h:48
~TKDTreeBinning()
Class&#39;s destructor.
void Initialize(unsigned int newPoints, unsigned int dim=1, ErrorType err=kValueError)
Definition: BinData.cxx:345
void Sort(Index n, const Element *a, Index *index, Bool_t down=kTRUE)
Definition: TMath.h:1151
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:9
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:33
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:873
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:53
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:104
const Bool_t kFALSE
Definition: RtypesCore.h:92
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:336
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:418
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.
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: RtypesCore.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:859
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.