Logo ROOT   6.10/09
Reference Guide
testkdTreeBinning.cxx
Go to the documentation of this file.
1 // ------------------------------------------------------------------------
2 //
3 // kdTreeBinning tutorial: bin the data in cells of equal content using a kd-tree
4 //
5 // Using TKDTree wrapper class as a data binning structure
6 // Plot the 2D data using the TH2Poly class
7 //
8 //
9 // Author: Bartolomeu Rabacal 11/2010
10 //
11 // ------------------------------------------------------------------------
12 
13 #include <cmath>
14 
15 #include "TKDTreeBinning.h"
16 #include "TH2D.h"
17 #include "TH2Poly.h"
18 #include "TStyle.h"
19 #include "TRandom3.h"
20 #include "TCanvas.h"
21 #include "TApplication.h"
22 #include "TMarker.h"
23 #include <iostream>
24 
25 bool verbose = false;
26 bool showGraphics = false;
27 
28 using std::cout;
29 using std::cerr;
30 using std::endl;
31 
33 
34  // -----------------------------------------------------------------------------------------------
35  // C r e a t e r a n d o m s a m p l e
36  // -----------------------------------------------------------------------------------------------
37 
38  const UInt_t DATASZ = 10000;
39  const UInt_t DATADIM = 2;
40  const UInt_t NBINS = 50;
41 
42  Double_t smp[DATASZ * DATADIM];
43 
44  double mu[2] = {0,2};
45  double sig[2] = {2,3};
46  TRandom3 r;
47  r.SetSeed(1);
48  for (UInt_t i = 0; i < DATADIM; ++i)
49  for (UInt_t j = 0; j < DATASZ; ++j)
50  smp[DATASZ * i + j] = r.Gaus(mu[i], sig[i]);
51 
52  UInt_t h1bins = (UInt_t) std::sqrt(double(NBINS));
53 
54  TH2D* h1 = new TH2D("h1BinTest", "Regular binning", h1bins, -5., 5., h1bins, -5., 5.);
55  for (UInt_t j = 0; j < DATASZ; ++j)
56  h1->Fill(smp[j], smp[DATASZ + j]);
57 
58 
59  // ---------------------------------------------------------------------------------------------
60  // C r e a t e K D T r e e B i n n i n g o b j e c t w i t h T H 2 P o l y p l o t t i n g
61  // ---------------------------------------------------------------------------------------------
62 
63  TKDTreeBinning* kdBins = new TKDTreeBinning(DATASZ, DATADIM, smp, NBINS);
64 
65  UInt_t nbins = kdBins->GetNBins();
66  UInt_t dim = kdBins->GetDim();
67 
68  const Double_t* binsMinEdges = kdBins->GetBinsMinEdges();
69  const Double_t* binsMaxEdges = kdBins->GetBinsMaxEdges();
70 
71  TH2Poly* h2pol = new TH2Poly("h2PolyBinTest", "KDTree binning", kdBins->GetDataMin(0), kdBins->GetDataMax(0), kdBins->GetDataMin(1), kdBins->GetDataMax(1));
72 
73  for (UInt_t i = 0; i < nbins; ++i) {
74  UInt_t edgeDim = i * dim;
75  h2pol->AddBin(binsMinEdges[edgeDim], binsMinEdges[edgeDim + 1], binsMaxEdges[edgeDim], binsMaxEdges[edgeDim + 1]);
76  }
77 
78  for (UInt_t i = 1; i <= kdBins->GetNBins(); ++i)
79  h2pol->SetBinContent(i, kdBins->GetBinDensity(i - 1));
80 
81  int ibinMin = kdBins->GetBinMinDensity();
82  int ibinMax = kdBins->GetBinMaxDensity();
83 
84  std::cout << "Bin with minimum density: " << ibinMin << " density = " << kdBins->GetBinDensity(ibinMin) << " content = " << kdBins->GetBinContent(ibinMin) << std::endl;
85  std::cout << "Bin with maximum density: " << ibinMax << " density = " << kdBins->GetBinDensity(ibinMax) << " content = " << kdBins->GetBinContent(ibinMin) << std::endl;
86 
87  if (kdBins->GetBinDensity(ibinMax) != DATASZ/NBINS)
88  Error("testkdTreeBinning","Wrong bin content");
89 
90  // order bins by density
91  kdBins->SortBinsByDensity(true);
92 
93  if (kdBins->GetBinMinDensity() != 0) Error("testkdTreeBinning","Wrong minimum bin after sorting");
94  if (kdBins->GetBinMaxDensity() != nbins-1) Error("testkdTreeBinning","Wrong maximum bin after sorting");
95 
96  if (showGraphics) {
97  new TCanvas();
98  h2pol->Draw("COLZ L");
99  gPad->Update();
100  }
101 
102 
103  // test find bin
104  int ntimes = (verbose) ? 2 : 200;
105  double point[2] = {0,0};
106 // double binCenter[2];
107  gRandom->SetSeed(0);
108  bool ok = true;
109  for (int itimes = 0; itimes < ntimes; itimes++) {
110 
111  // generate a random point in 2D
112  point[0] = gRandom->Uniform(-5,5);
113  point[1] = gRandom->Uniform(-5,5);
114  // int inode = tree->FindNode(point);
115  // inode = inode - tree->GetNNodes();
116 
117  int ibin = kdBins->FindBin(point);
118 
119  const double * binCenter = kdBins->GetBinCenter(ibin);
120  const double * binMin = kdBins->GetBinMinEdges(ibin);
121  const double * binMax = kdBins->GetBinMaxEdges(ibin);
122  if (binCenter) {
123 
124  if (verbose) {
125  std::cout << "**** point *** " << itimes << "\n";
126  std::cout << " point x " << point[0] << " BIN CENTER is " << binCenter[0] << " min " << binMin[0] << " max " << binMax[0] << std::endl;
127  std::cout << " point y " << point[1] << " BIN CENTER is " << binCenter[1] << " min " << binMin[1] << " max " << binMax[1] << std::endl;
128  }
129 
130  ok &= point[0] > binMin[0] && point[0] < binMax[0];
131  ok &= point[1] > binMin[1] && point[1] < binMax[1];
132  if (!ok) {
133  Error ("testkdTreeBinning::FindBin"," Point is not in the right bin " );
134  std::cout << " point x " << point[0] << " BIN CENTER is " << binCenter[0] << " min " << binMin[0] << " max " << binMax[0] << std::endl;
135  std::cout << " point y " << point[1] << " BIN CENTER is " << binCenter[1] << " min " << binMin[1] << " max " << binMax[1] << std::endl;
136  }
137 
138  if (itimes < 2 && showGraphics ) {
139 
140  TMarker * m1 = new TMarker(point[0], point[1],20 );
141  TMarker * m2 = new TMarker(binCenter[0], binCenter[1], 21);
142  m1->Draw();
143  m2->Draw();
144  }
145  }
146  else
147  Error("testkdTreeBinning::FindBin"," Bin %d is not existing ",ibin);
148 
149 
150  }
151 
152  return;
153 
154 }
155 
156 int main(int argc, char **argv)
157 {
158  // Parse command line arguments
159  for (Int_t i=1 ; i<argc ; i++) {
160  std::string arg = argv[i] ;
161  if (arg == "-g") {
162  showGraphics = true;
163  }
164  if (arg == "-v") {
165  showGraphics = true;
166  verbose = true;
167  }
168  if (arg == "-h") {
169  cerr << "Usage: " << argv[0] << " [-g] [-v]\n";
170  cerr << " where:\n";
171  cerr << " -g : graphics mode\n";
172  cerr << " -v : verbose mode";
173  cerr << endl;
174  return -1;
175  }
176  }
177 
178  TApplication* theApp = 0;
179  if ( showGraphics )
180  theApp = new TApplication("App",&argc,argv);
181 
183 
184  if ( showGraphics )
185  {
186  theApp->Run();
187  delete theApp;
188  theApp = 0;
189  }
190 
191  return 0;
192 
193 }
194 
Random number generator class based on M.
Definition: TRandom3.h:27
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}.
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...
int main(int argc, char **argv)
virtual Double_t Gaus(Double_t mean=0, Double_t sigma=1)
Samples a random number from the standard Normal (Gaussian) Distribution with the given mean and sigm...
Definition: TRandom.cxx:235
UInt_t GetBinMinDensity() const
Return the bin with minimum density.
bool verbose
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}.
virtual void Draw(Option_t *option="")
Draw this marker with its current attributes.
Definition: TMarker.cxx:143
virtual void SetSeed(ULong_t seed=0)
Set the random generator sequence if seed is 0 (default value) a TUUID is generated and used to fill ...
Definition: TRandom3.cxx:201
UInt_t GetNBins() const
Returns the number of bins.
Manages Markers.
Definition: TMarker.h:23
int Int_t
Definition: RtypesCore.h:41
int nbins[3]
double sqrt(double)
virtual void Run(Bool_t retrn=kFALSE)
Main application eventloop. Calls system dependent eventloop via gSystem.
virtual void SetSeed(ULong_t seed=0)
Set the random generator seed.
Definition: TRandom.cxx:568
virtual Int_t AddBin(TObject *poly)
Adds a new bin to the histogram.
Definition: TH2Poly.cxx:214
<- TKDTreeBinning - A class providing multidimensional binning ->
TH1F * h1
Definition: legend1.C:5
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...
void Error(const char *location, const char *msgfmt,...)
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.
Double_t GetBinDensity(UInt_t bin) const
Returns the density in bin.
TRandom2 r(17)
virtual void Draw(Option_t *option="")
Draw this histogram with options.
Definition: TH1.cxx:2851
UInt_t GetBinContent(UInt_t bin) const
Returns the number of points in bin. &#39;bin&#39; is between 0 and fNBins - 1.
void SetBinContent(Int_t bin, Double_t content)
Sets the contents of the input bin to the input content Negative values between -1 and -9 are for the...
Definition: TH2Poly.cxx:1247
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
void testkdTreeBinning()
R__EXTERN TRandom * gRandom
Definition: TRandom.h:62
void SortBinsByDensity(Bool_t sortAsc=kTRUE)
Sorts bins by their density.
The Canvas class.
Definition: TCanvas.h:31
double Double_t
Definition: RtypesCore.h:55
UInt_t GetBinMaxDensity() const
Return the bin with maximum density.
THist< 2, double, THistStatContent, THistStatUncertainty > TH2D
Definition: THist.hxx:316
virtual Double_t Uniform(Double_t x1=1)
Returns a uniform deviate on the interval (0, x1).
Definition: TRandom.cxx:606
const Double_t * GetBinMaxEdges(UInt_t bin) const
Returns the bin&#39;s maximum edges. &#39;bin&#39; is between 0 and fNBins - 1.
#define gPad
Definition: TVirtualPad.h:284
bool showGraphics
This class creates the ROOT Application Environment that interfaces to the windowing system eventloop...
Definition: TApplication.h:39
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.
Int_t Fill(Double_t)
Invalid Fill method.
Definition: TH2.cxx:292
2D Histogram with Polygonal Bins
Definition: TH2Poly.h:66
tomato 2-D histogram with a double per channel (see TH1 documentation)}
Definition: TH2.h:290