Logo ROOT  
Reference Guide
 
Loading...
Searching...
No Matches
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 "TBuffer.h"
22
24
25/**
26\class TKDTreeBinning
27<- TKDTreeBinning - A class providing multidimensional binning ->
28
29The class implements multidimensional binning by constructing a TKDTree inner structure from the
30data which is used as the bins.
31The bins are retrieved as two double*, one for the minimum bin edges,
32the other as the maximum bin edges. For one dimension one of these is enough to correctly define the bins.
33For the multidimensional case both minimum and maximum ones are necessary for the bins to be well defined.
34The bin edges of d-dimensional data is a d-tet of the bin's thresholds. For example if d=3 the minimum bin
35edges of bin b is of the form of the following array: {xbmin, ybmin, zbmin}.
36You also have the possibility to sort the bins by their density.
37
38Details of usage can be found in `$ROOTSYS/tutorials/math/kdTreeBinning.C` and more information on
39the embedded TKDTree documentation.
40
41@ingroup MathCore
42
43*/
44
46 // Boolean functor whose predicate depends on the bin's density. Used for ascending sort.
47 CompareAsc(const TKDTreeBinning* treebins) : bins(treebins) {}
49 return bins->GetBinDensity(bin1) < bins->GetBinDensity(bin2);
50 }
52};
53
55 // Boolean functor whose predicate depends on the bin's density. Used for descending sort.
56 CompareDesc(const TKDTreeBinning* treebins) : bins(treebins) {}
58 return bins->GetBinDensity(bin1) > bins->GetBinDensity(bin2);
59 }
61};
62
63/// Class's constructor taking the size of the data points, dimension, a data array and the number
64/// of bins (default = 100). It is reccomended to have the number of bins as an exact divider of
65/// the data size.
66/// The data array must be organized with a stride=1 for the points and = N (the dataSize) for the dimension.
67///
68/// Thus data[] = x1,x2,x3,......xN, y1,y2,y3......yN, z1,z2,...........zN,....
69///
70/// Note that the passed dataSize is not the size of the array but is the number of points (N)
71/// The size of the array must be at least dataDim*dataSize
72///
73TKDTreeBinning::TKDTreeBinning(UInt_t dataSize, UInt_t dataDim, Double_t* data, UInt_t nBins, bool adjustBinEdges)
74: fData(0), fBinMinEdges(std::vector<Double_t>()), fBinMaxEdges(std::vector<Double_t>()), fDataBins((TKDTreeID*)0), fDim(dataDim),
75fDataSize(dataSize), fDataThresholds(std::vector<std::pair<Double_t, Double_t> >(fDim, std::make_pair(0., 0.))),
76fIsSorted(kFALSE), fIsSortedAsc(kFALSE), fBinsContent(std::vector<UInt_t>()) {
77 if (adjustBinEdges) SetBit(kAdjustBinEdges);
78 if (data) {
79 SetData(data);
80 SetNBins(nBins);
81 } else {
82 if (fData.empty())
83 this->Warning("TKDTreeBinning", "Data is nil. Nothing is built.");
84 }
85}
86/// Class's constructor taking the size of the data points, dimension, a data vector and the number
87/// of bins (default = 100). It is reccomended to have the number of bins as an exact divider of
88/// the data size.
89/// The data array must be organized with a stride=1 for the points and = N (the dataSize) for the dimension.
90///
91/// Thus data[] = x1,x2,x3,......xN, y1,y2,y3......yN, z1,z2,...........zN,....
92///
93/// Note that the passed data vector may contains a larger size, in case extra coordinates are associated but not used
94/// in building the kdtree
95/// The size of thedata vector must be at least dataDim*dataSize
96///
97TKDTreeBinning::TKDTreeBinning(UInt_t dataSize, UInt_t dataDim, const std::vector<double> &data, UInt_t nBins, bool adjustBinEdges)
98: fData(0), fBinMinEdges(std::vector<Double_t>()), fBinMaxEdges(std::vector<Double_t>()), fDataBins((TKDTreeID*)0), fNBins (nBins), fDim(dataDim),
99fDataSize(dataSize), fDataThresholds(std::vector<std::pair<Double_t, Double_t> >(fDim, std::make_pair(0., 0.))),
100fIsSorted(kFALSE), fIsSortedAsc(kFALSE), fBinsContent(std::vector<UInt_t>()) {
101 if (adjustBinEdges) SetBit(kAdjustBinEdges);
102 if (!data.empty()) {
103 SetData(data);
104 SetNBins(nBins);
105 } else {
106 if (fData.empty())
107 this->Warning("TKDTreeBinning", "Data is nil. Nothing is built.");
108 }
109}
110
111/// Default constructor (for I/O)
113// fData(nullptr),
114 fDataBins(nullptr),
115 fNBins (0),
116 fDim(0),
117 fDataSize(0),
118 fIsSorted(kFALSE), fIsSortedAsc(kFALSE)
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();
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
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}
223void 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)
239}
240
242 // Sets the bins' content
243 fBinsContent.resize(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());
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
369std::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
405std::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
504const 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
540
541 // Add also the upper(lower) edge to the min (max) list
542 if (sortAsc) {
543 fBinMinEdges.push_back(fBinMaxEdges.back());
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,
553/// Array must be deleted as "delete [] res"
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 nullptr;
566}
567
568/// Returns a pointer to the vector of the bin widths. 'bin' is between 0 and fNBins - 1
569/// Array must be deleted as "delete [] res"
571 if(bin < fNBins) {
572 Double_t* result = new Double_t[fDim];
573 std::pair<const Double_t*, const Double_t*> binEdges = GetBinEdges(bin);
574 for (UInt_t i = 0; i < fDim; ++i) {
575 result[i] = (binEdges.second[i] - binEdges.first[i]);
576 }
577 return result;
578 }
579 this->Warning("GetBinWidth", "No such bin. Returning null pointer.");
580 this->Info("GetBinWidth", "'bin' is between 0 and %d.", fNBins - 1);
581 return nullptr;
582}
583
584/// Return the bin with maximum density
586 if (fIsSorted) {
587 if (fIsSortedAsc)
588 return fNBins - 1;
589 else return 0;
590 }
591 UInt_t* indices = new UInt_t[fNBins];
592 for (UInt_t i = 0; i < fNBins; ++i)
593 indices[i] = i;
594 UInt_t result = *std::max_element(indices, indices + fNBins, CompareAsc(this));
595 delete [] indices;
596 return result;
597}
598
599/// Return the bin with minimum density
601 if (fIsSorted) {
602 if (!fIsSortedAsc)
603 return fNBins - 1;
604 else return 0;
605 }
606 UInt_t* indices = new UInt_t[fNBins];
607 for (UInt_t i = 0; i < fNBins; ++i)
608 indices[i] = i;
609 UInt_t result = *std::min_element(indices, indices + fNBins, CompareAsc(this));
610 delete [] indices;
611 return result;
612}
613
614/// Fill the bin data set (class ROOT::Fit::BinData) with the result of the TKDTree binning
616 if (!fDataBins) return;
617 data.Initialize(fNBins, fDim);
618 for (unsigned int i = 0; i < fNBins; ++i) {
619 data.Add( GetBinMinEdges(i), GetBinDensity(i), std::sqrt(double(GetBinContent(i) ))/ GetBinVolume(i) );
621 }
622}
623
624/// find the corresponding bin index given the coordinate of a point
626
627 Int_t inode = fDataBins->FindNode(point);
628 // find node return the index in the total nodes and the bins are only the terminal ones
629 // so we subtract all the non-terminal nodes
630 inode -= fDataBins->GetNNodes();
631 R__ASSERT( inode >= 0);
632 UInt_t bin = inode;
633
634 if (!fIsSorted) return bin;
635 //return std::distance(fIndices.begin(), std::find(fIndices.begin(), fIndices.end(), bin ) );
636 return fIndices[bin];
637}
638
639/// Return the corresponding point belonging to the bin i
640std::vector<std::vector<Double_t> > TKDTreeBinning::GetPointsInBin(UInt_t bin) const {
641 std::vector<Double_t> point(fDim);
642 std::vector< std::vector<Double_t> > thePoints;
643 if (fData.size() == 0) {
644 Error("GetPointsInBin","Internal data set is not valid");
645 return thePoints;
646 }
647 if (!fDataBins) {
648 Error("GetPointsInBin","Internal TKDTree is not valid");
649 return thePoints;
650 }
651 if (bin >= fNBins) {
652 Error("GetPointsInBin","Invalid bin number");
653 return thePoints;
654 }
655 // need to find the bin number including the non-terminal node
656 int inode = bin + fDataBins->GetNNodes();
657 auto indices = fDataBins->GetPointsIndexes(inode);
658 int npoints = fDataBins->GetNPointsNode(inode);
659 thePoints.resize(npoints);
660 for (int ipoint = 0; ipoint < npoints; ++ipoint) {
661 for (unsigned int idim = 0; idim < fDim; ++idim) {
662 point[idim] = fData[idim*fDataSize+indices[ipoint] ];
663 }
664 thePoints[ipoint] = point;
665 }
666 return thePoints;
667}
668
669////////////////////////////////////////////////////////////////////////////////
670/// Stream a class object.
671void TKDTreeBinning::Streamer(TBuffer &b) {
672 if (b.IsReading() ) {
673 UInt_t R__s, R__c;
674 Version_t v = b.ReadVersion(&R__s, &R__c);
675 b.ReadClassBuffer(TKDTreeBinning::Class(), this, v, R__s, R__c);
676 // we need to rebuild the TKDTree
677 if (fDataBins) delete fDataBins;
679 }
680 else {
681 // case of writing
682 b.WriteClassBuffer(TKDTreeBinning::Class(),this);
683 //std::cout << "writing npar = " << GetNpar() << std::endl;
684 }
685}
#define b(i)
Definition RSha256.hxx:100
short Version_t
Definition RtypesCore.h:65
unsigned int UInt_t
Definition RtypesCore.h:46
const Bool_t kFALSE
Definition RtypesCore.h:101
double Double_t
Definition RtypesCore.h:59
const Bool_t kTRUE
Definition RtypesCore.h:100
#define ClassImp(name)
Definition Rtypes.h:364
#define R__ASSERT(e)
Definition TError.h:118
TKDTree< Int_t, Double_t > TKDTreeID
Definition TKDTree.h:104
Class describing the binned data sets : vectors of x coordinates, y values and optionally error on y ...
Definition BinData.h:52
void AddBinUpEdge(const double *xup)
add the bin width data, a pointer to an array with the bin upper edge information.
Definition BinData.cxx:621
void Add(double x, double y)
add one dim data with only coordinate and values
Definition BinData.cxx:416
void Initialize(unsigned int newPoints, unsigned int dim=1, ErrorType err=kValueError)
Definition BinData.cxx:343
Buffer base class used for serializing objects.
Definition TBuffer.h:43
<- TKDTreeBinning - A class providing multidimensional binning ->
std::vector< UInt_t > fBinsContent
Flags if the bin edges are sorted densitywise (or by bin-edge for 1D) in ascending order.
void ReadjustMinBinEdges(Double_t *binEdges)
TKDTreeBinning()
Default constructor (for I/O)
std::pair< const Double_t *, const Double_t * > GetBinsEdges() const
Returns a pair of an array with all bins minimum and maximum edges.
const Double_t * GetBinMaxEdges(UInt_t bin) const
Returns the bin's maximum edges. 'bin' is between 0 and fNBins - 1.
const Double_t * GetBinWidth(UInt_t bin) const
Returns a pointer to the vector of the bin widths.
const Double_t * GetBinsMaxEdges() const
Returns an array with all bins' maximum edges The edges are arranges as xmax_1,ymax_1,...
void SetCommonBinEdges(Double_t *binEdges)
const Double_t * GetBinCenter(UInt_t bin) const
Returns the geometric center of of the bin.
UInt_t GetBinMaxDensity() const
Return the bin with maximum density.
~TKDTreeBinning()
Class's destructor.
std::pair< const Double_t *, const Double_t * > GetBinEdges(UInt_t bin) const
Returns a pir with the bin's edges. 'bin' is between 0 and fNBins - 1.
UInt_t GetNBins() const
Returns the number of bins.
void SetData(Double_t *data)
Disallowed assign operator.
void SetNBins(UInt_t bins)
Sets binning inner structure.
UInt_t fDataSize
The data dimension.
Double_t GetDataMax(UInt_t dim) const
Returns the maximum of the data in the dim coordinate. 'dim' is between 0 and fDim - 1.
const Double_t * GetBinMinEdges(UInt_t bin) const
Returns the bin's minimum edges. 'bin' is between 0 and fNBins - 1.
Double_t GetDataMin(UInt_t dim) const
Returns the minimum of the data in the dim coordinate. 'dim' is between 0 and fDim - 1.
UInt_t GetBinContent(UInt_t bin) const
Returns the number of points in bin. 'bin' is between 0 and fNBins - 1.
void FillBinData(ROOT::Fit::BinData &data) const
Fill the bin data set (class ROOT::Fit::BinData) with the result of the TKDTree binning.
std::vector< Double_t > fData
Double_t GetBinVolume(UInt_t bin) const
Returns the (hyper)volume of bin. 'bin' is between 0 and fNBins - 1.
std::vector< std::pair< Double_t, Double_t > > fDataThresholds
The data size.
void SortBinsByDensity(Bool_t sortAsc=kTRUE)
Sorts bins by their density.
std::vector< UInt_t > fIndices
The maximum values for the bins' edges for each dimension.
std::vector< Double_t > fBinMinEdges
[fDataSize*fDim] The data from which a KDTree partition is computed for binning
std::vector< Double_t > fBinMaxEdges
The minimum values for the bins' edges for each dimension.
void ReadjustMaxBinEdges(Double_t *binEdges)
std::vector< std::vector< std::pair< Bool_t, Bool_t > > > fCheckedBinEdges
Minimum and maximum data values.
UInt_t fDim
The number of bins.
UInt_t FindBin(const Double_t *point) const
find the corresponding bin index given the coordinate of a point
Double_t GetBinDensity(UInt_t bin) const
Returns the density in bin.
TKDTreeID * fDataBins
Index of the bins in the kd-tree (needed when bins are sorted)
Bool_t fIsSortedAsc
Flags if the bin edges are sorted densitywise (or by bin endges in case of 1-dim )
std::vector< std::vector< Double_t > > GetPointsInBin(UInt_t bin) const
Return the corresponding point belonging to the bin i.
std::vector< std::map< Double_t, std::vector< UInt_t > > > fCommonBinEdges
! Auxiliary structure for readjusting the bin edges. Keeps the common bin boundaries
const Double_t * GetBinsMinEdges() const
Returns an array with all bins' minimum edges The edges are arranges as xmin_1,ymin_1,...
void SetBinMinMaxEdges(Double_t *binEdges)
UInt_t GetBinMinDensity() const
Return the bin with minimum density.
const Double_t * GetOneDimBinEdges() const
Returns a pointer to the vector of the bin edges for one dimensional binning only.
const Double_t * GetDimData(UInt_t dim) const
UInt_t GetDim() const
Returns the number of dimensions.
TKDTreeID * GetTree() const
Returns the kD-Tree structure of the binning.
const Double_t * SortOneDimBinEdges(Bool_t sortAsc=kTRUE)
Sort the one-dimensional bin edges and retuns a pointer to them.
Class implementing a kd-tree.
Definition TKDTree.h:10
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
Index GetBucketSize()
Definition TKDTree.h:48
Index FindNode(const Value *point) const
returns the index of the terminal node to which point belongs (index in the fAxis,...
Definition TKDTree.cxx:677
Index * GetPointsIndexes(Int_t node) const
return the indices of the points in that terminal node for all the nodes except last,...
Definition TKDTree.cxx:817
void Build()
Build the kd-tree.
Definition TKDTree.cxx:411
Int_t GetNNodes() const
Definition TKDTree.h:33
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
Value * GetBoundary(const Int_t node)
Get a boundary.
Definition TKDTree.cxx:1211
R__ALWAYS_INLINE Bool_t TestBit(UInt_t f) const
Definition TObject.h:201
virtual void Warning(const char *method, const char *msgfmt,...) const
Issue warning message.
Definition TObject.cxx:949
void SetBit(UInt_t f, Bool_t set)
Set or unset the user status bits as specified in f.
Definition TObject.cxx:766
virtual void Error(const char *method, const char *msgfmt,...) const
Issue error message.
Definition TObject.cxx:963
virtual void Info(const char *method, const char *msgfmt,...) const
Issue info message.
Definition TObject.cxx:937
void Sort(Index n, const Element *a, Index *index, Bool_t down=kTRUE)
Definition TMathBase.h:358
Definition first.py:1
const TKDTreeBinning * bins
CompareAsc(const TKDTreeBinning *treebins)
Bool_t operator()(UInt_t bin1, UInt_t bin2)
const TKDTreeBinning * bins
Bool_t operator()(UInt_t bin1, UInt_t bin2)
CompareDesc(const TKDTreeBinning *treebins)