Logo ROOT  
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
30The class implements multidimensional binning by constructing a TKDTree inner structure from the
31data which is used as the bins.
32The bins are retrieved as two double*, one for the minimum bin edges,
33the other as the maximum bin edges. For one dimension one of these is enough to correctly define the bins.
34For the multidimensional case both minimum and maximum ones are necessary for the bins to be well defined.
35The bin edges of d-dimensional data is a d-tet of the bin's thresholds. For example if d=3 the minimum bin
36edges of bin b is of the form of the following array: {xbmin, ybmin, zbmin}.
37You also have the possibility to sort the bins by their density.
38
39Details of usage can be found in `$ROOTSYS/tutorials/math/kdTreeBinning.C` and more information on
40the 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///
74TKDTreeBinning::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),
76fDataSize(dataSize), fDataThresholds(std::vector<std::pair<Double_t, Double_t> >(fDim, std::make_pair(0., 0.))),
77fIsSorted(kFALSE), fIsSortedAsc(kFALSE), fBinsContent(std::vector<UInt_t>()) {
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///
98TKDTreeBinning::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),
100fDataSize(dataSize), fDataThresholds(std::vector<std::pair<Double_t, Double_t> >(fDim, std::make_pair(0., 0.))),
101fIsSorted(kFALSE), fIsSortedAsc(kFALSE), fBinsContent(std::vector<UInt_t>()) {
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),
119 fIsSorted(kFALSE), fIsSortedAsc(kFALSE)
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();
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
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}
224void 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)
240}
241
243 // Sets the bins' content
244 fBinsContent.resize(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());
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;
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];
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
370std::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
406std::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
505const 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
541
542 // Add also the upper(lower) edge to the min (max) list
543 if (sortAsc) {
544 fBinMinEdges.push_back(fBinMaxEdges.back());
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,
554/// Array must be deleted as "delete [] res"
556 if(bin < fNBins) {
557 Double_t* result = new Double_t[fDim];
558 std::pair<const Double_t*, const Double_t*> binEdges = GetBinEdges(bin);
559 for (UInt_t i = 0; i < fDim; ++i) {
560 result[i] = (binEdges.second[i] + binEdges.first[i]) / 2.;
561 }
562 return result;
563 }
564 this->Warning("GetBinCenter", "No such bin. Returning null pointer.");
565 this->Info("GetBinCenter", "'bin' is between 0 and %d.", fNBins - 1);
566 return nullptr;
567}
568
569/// Returns a pointer to the vector of the bin widths. 'bin' is between 0 and fNBins - 1
570/// Array must be deleted as "delete [] res"
572 if(bin < fNBins) {
573 Double_t* result = new Double_t[fDim];
574 std::pair<const Double_t*, const Double_t*> binEdges = GetBinEdges(bin);
575 for (UInt_t i = 0; i < fDim; ++i) {
576 result[i] = (binEdges.second[i] - binEdges.first[i]);
577 }
578 return result;
579 }
580 this->Warning("GetBinWidth", "No such bin. Returning null pointer.");
581 this->Info("GetBinWidth", "'bin' is between 0 and %d.", fNBins - 1);
582 return nullptr;
583}
584
585/// Return the bin with maximum density
587 if (fIsSorted) {
588 if (fIsSortedAsc)
589 return fNBins - 1;
590 else return 0;
591 }
592 UInt_t* indices = new UInt_t[fNBins];
593 for (UInt_t i = 0; i < fNBins; ++i)
594 indices[i] = i;
595 UInt_t result = *std::max_element(indices, indices + fNBins, CompareAsc(this));
596 delete [] indices;
597 return result;
598}
599
600/// Return the bin with minimum density
602 if (fIsSorted) {
603 if (!fIsSortedAsc)
604 return fNBins - 1;
605 else return 0;
606 }
607 UInt_t* indices = new UInt_t[fNBins];
608 for (UInt_t i = 0; i < fNBins; ++i)
609 indices[i] = i;
610 UInt_t result = *std::min_element(indices, indices + fNBins, CompareAsc(this));
611 delete [] indices;
612 return result;
613}
614
615/// Fill the bin data set (class ROOT::Fit::BinData) with the result of the TKDTree binning
617 if (!fDataBins) return;
618 data.Initialize(fNBins, fDim);
619 for (unsigned int i = 0; i < fNBins; ++i) {
620 data.Add( GetBinMinEdges(i), GetBinDensity(i), std::sqrt(double(GetBinContent(i) ))/ GetBinVolume(i) );
622 }
623}
624
625/// find the corresponding bin index given the coordinate of a point
627
628 Int_t inode = fDataBins->FindNode(point);
629 // find node return the index in the total nodes and the bins are only the terminal ones
630 // so we subtract all the non-terminal nodes
631 inode -= fDataBins->GetNNodes();
632 R__ASSERT( inode >= 0);
633 UInt_t bin = inode;
634
635 if (!fIsSorted) return bin;
636 //return std::distance(fIndices.begin(), std::find(fIndices.begin(), fIndices.end(), bin ) );
637 return fIndices[bin];
638}
639
640/// Return the corresponding point belonging to the bin i
641std::vector<std::vector<Double_t> > TKDTreeBinning::GetPointsInBin(UInt_t bin) const {
642 std::vector<Double_t> point(fDim);
643 std::vector< std::vector<Double_t> > thePoints;
644 if (fData.size() == 0) {
645 Error("GetPointsInBin","Internal data set is not valid");
646 return thePoints;
647 }
648 if (!fDataBins) {
649 Error("GetPointsInBin","Internal TKDTree is not valid");
650 return thePoints;
651 }
652 if (bin >= fNBins) {
653 Error("GetPointsInBin","Invalid bin number");
654 return thePoints;
655 }
656 // need to find the bin number including the non-terminal node
657 int inode = bin + fDataBins->GetNNodes();
658 auto indices = fDataBins->GetPointsIndexes(inode);
659 int npoints = fDataBins->GetNPointsNode(inode);
660 thePoints.resize(npoints);
661 for (int ipoint = 0; ipoint < npoints; ++ipoint) {
662 for (unsigned int idim = 0; idim < fDim; ++idim) {
663 point[idim] = fData[idim*fDataSize+indices[ipoint] ];
664 }
665 thePoints[ipoint] = point;
666 }
667 return thePoints;
668}
669
670////////////////////////////////////////////////////////////////////////////////
671/// Stream a class object.
672void TKDTreeBinning::Streamer(TBuffer &b) {
673 if (b.IsReading() ) {
674 UInt_t R__s, R__c;
675 Version_t v = b.ReadVersion(&R__s, &R__c);
676 b.ReadClassBuffer(TKDTreeBinning::Class(), this, v, R__s, R__c);
677 // we need to rebuild the TKDTree
678 if (fDataBins) delete fDataBins;
680 }
681 else {
682 // case of writing
683 b.WriteClassBuffer(TKDTreeBinning::Class(),this);
684 //std::cout << "writing npar = " << GetNpar() << std::endl;
685 }
686}
void Class()
Definition: Class.C:29
#define b(i)
Definition: RSha256.hxx:100
short Version_t
Definition: RtypesCore.h:63
unsigned int UInt_t
Definition: RtypesCore.h:44
const Bool_t kFALSE
Definition: RtypesCore.h:90
bool Bool_t
Definition: RtypesCore.h:61
double Double_t
Definition: RtypesCore.h:57
const Bool_t kTRUE
Definition: RtypesCore.h:89
#define ClassImp(name)
Definition: Rtypes.h:361
#define R__ASSERT(e)
Definition: TError.h:96
TKDTree< Int_t, Double_t > TKDTreeID
Definition: TKDTree.h:104
double sqrt(double)
TRObject operator()(const T1 &t1) const
Class describing the binned data sets : vectors of x coordinates, y values and optionally error on y ...
Definition: BinData.h:53
void AddBinUpEdge(const double *xup)
add the bin width data, a pointer to an array with the bin upper edge information.
Definition: BinData.cxx:627
void Add(double x, double y)
add one dim data with only coordinate and values
Definition: BinData.cxx:422
void Initialize(unsigned int newPoints, unsigned int dim=1, ErrorType err=kValueError)
Definition: BinData.cxx:349
Buffer base class used for serializing objects.
Definition: TBuffer.h:42
<- 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.
friend struct CompareDesc
! Predicate for descending sort
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.
friend struct CompareAsc
! Predicate for ascending sort
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:187
virtual void Warning(const char *method, const char *msgfmt,...) const
Issue warning message.
Definition: TObject.cxx:877
void SetBit(UInt_t f, Bool_t set)
Set or unset the user status bits as specified in f.
Definition: TObject.cxx:694
virtual void Error(const char *method, const char *msgfmt,...) const
Issue error message.
Definition: TObject.cxx:891
virtual void Info(const char *method, const char *msgfmt,...) const
Issue info message.
Definition: TObject.cxx:865
VecExpr< UnaryOp< Fabs< T >, VecExpr< A, T, D >, T >, T, D > fabs(const VecExpr< A, T, D > &rhs)
static constexpr double second
void Sort(Index n, const Element *a, Index *index, Bool_t down=kTRUE)
Definition: TMathBase.h:362
Definition: first.py:1
REAL epsilon
Definition: triangle.c:617