Logo ROOT  
Reference Guide
 
Loading...
Searching...
No Matches
SparseData.cxx
Go to the documentation of this file.
1// @(#)root/mathcore:$Id$
2// Author: David Gonzalez Maline Wed Aug 28 15:33:03 2009
3
4/**********************************************************************
5 * *
6 * Copyright (c) 2006 LCG ROOT Math Team, CERN/PH-SFT *
7 * *
8 * *
9 **********************************************************************/
10
11// Implementation file for class BinData
12
13#include <iostream>
14#include <iterator>
15#include <algorithm>
16
17#include <vector>
18#include <list>
19
20#include <stdexcept>
21
22#include <cmath>
23#include <limits>
24
25// #include "TMath.h"
26#include "Fit/BinData.h"
27#include "Fit/SparseData.h"
28
29using namespace std;
30
31namespace ROOT {
32
33 namespace Fit {
34
35 //This class is a helper. It represents a bin in N
36 //dimensions. The change in the name is to avoid name collision.
37 class Box
38 {
39 public:
40 // Creates a Box with limits specified by the vectors and
41 // content=value and error=error
42 Box(const vector<double>& min, const vector<double>& max,
43 const double value = 0.0, const double error = 1.0):
44 fMin(min), fMax(max), fVal(value), fError(error)
45 { }
46
47 // Compares to Boxes to see if they are equal in all its
48 // variables. This is to be used by the std::find algorithm
49 bool operator==(const Box& b)
50 { return (fMin == b.fMin) && (fMax == b.fMax)
51 && (fVal == b.fVal) && (fError == b.fError); }
52
53 // Get the list of minimum coordinates
54 const vector<double>& GetMin() const { return fMin; }
55 // Get the list of maximum coordinates
56 const vector<double>& GetMax() const { return fMax; }
57 // Get the value of the Box
58 double GetVal() const { return fVal; }
59 // Get the rror of the Box
60 double GetError() const { return fError; }
61
62 // Add an amount to the content of the Box
63 void AddVal(const double value) { fVal += value; }
64
65 friend class BoxContainer;
66 friend ostream& operator <<(ostream& os, const Box& b);
67
68 private:
69 vector<double> fMin;
70 vector<double> fMax;
71 double fVal;
72 double fError;
73 };
74
75 // This class is just a helper to be used in std::for_each to
76 // simplify the code later. It's just a definition of a method
77 // that will discern whether a Box is included into another one
79 {
80 private:
81 const Box& fBox;
82 public:
83 //Constructs the BoxContainer object with a Box that is meant
84 //to include another one that will be provided later
85 BoxContainer(const Box& b): fBox(b) {}
86
87 bool operator() (const Box& b1)
88 { return operator()(fBox, b1); }
89
90 // Looks if b2 is included in b1
91 bool operator() (const Box& b1, const Box& b2)
92 {
93 bool isIn = true;
94 vector<double>::const_iterator boxit = b2.fMin.begin();
95 vector<double>::const_iterator bigit = b1.fMax.begin();
96 while ( isIn && boxit != b2.fMin.end() )
97 {
98 if ( (*boxit) >= (*bigit) ) isIn = false;
99 ++boxit;
100 ++bigit;
101 }
102
103 boxit = b2.fMax.begin();
104 bigit = b1.fMin.begin();
105 while ( isIn && boxit != b2.fMax.end() )
106 {
107 if ( (*boxit) <= (*bigit) ) isIn = false;
108 ++boxit;
109 ++bigit;
110 }
111
112 return isIn;
113 }
114 };
115
116 // Another helper class to be used in std::for_each to simplify
117 // the code later. It implements the operator() to know if a
118 // specified Box is big enough to contain any 'space' inside.
120 {
121 public:
122 AreaComparer(vector<double>::iterator iter):
123 fThereIsArea(true),
124 fIter(iter),
125 fLimit(8 * std::numeric_limits<double>::epsilon())
126 {};
127
128 void operator() (double value)
129 {
130 if ( fabs(value- (*fIter)) < fLimit )
131// if ( TMath::AreEqualRel(value, (*fIter), fLimit) )
132 fThereIsArea = false;
133
134 ++fIter;
135 }
136
137 bool IsThereArea() { return fThereIsArea; }
138
139 private:
141 vector<double>::iterator fIter;
142 double fLimit;
143 };
144
145
146 // This is the key of the SparseData structure. This method
147 // will, by recursion, divide the area passed as an argument in
148 // min and max into pieces to insert the Box defined by bmin and
149 // bmax. It will do so from the highest dimension until it gets
150 // to 1 and create the corresponding boxes to divide the
151 // original space.
152 void DivideBox( const vector<double>& min, const vector<double>& max,
153 const vector<double>& bmin, const vector<double>& bmax,
154 const unsigned int size, const unsigned int n,
155 list<Box>& l, const double val, const double error)
156 {
157 vector<double> boxmin(min);
158 vector<double> boxmax(max);
159
160 boxmin[n] = min[n];
161 boxmax[n] = bmin[n];
162 if ( for_each(boxmin.begin(), boxmin.end(), AreaComparer(boxmax.begin())).IsThereArea() )
163 l.push_back(Box(boxmin, boxmax));
164
165 boxmin[n] = bmin[n];
166 boxmax[n] = bmax[n];
167 if ( n == 0 )
168 {
169 if ( for_each(boxmin.begin(), boxmin.end(), AreaComparer(boxmax.begin())).IsThereArea() )
170 l.push_back(Box(boxmin, boxmax, val, error));
171 }
172 else
173 DivideBox(boxmin, boxmax, bmin, bmax, size, n-1, l, val, error);
174
175 boxmin[n] = bmax[n];
176 boxmax[n] = max[n];
177 if ( for_each(boxmin.begin(), boxmin.end(), AreaComparer(boxmax.begin())).IsThereArea() )
178 l.push_back(Box(boxmin, boxmax));
179 }
180
182 {
183 public:
184 void PushBack(Box& box) { fProxy.push_back(box); }
185 list<Box>::iterator Begin() { return fProxy.begin(); }
186 list<Box>::iterator End() { return fProxy.end(); }
187 void Remove(list<Box>::iterator it) { fProxy.erase(it); }
188 list<Box>& GetList() { return fProxy; }
189 private:
190 list<Box> fProxy;
191 };
192
193
194 SparseData::SparseData(vector<double>& min, vector<double>& max)
195 {
196 // Creates a SparseData convering the range defined by min
197 // and max. For this it will create an empty Box for that
198 // range.
199 Box originalBox(min, max);
200 fList = new ProxyListBox();
201 fList->PushBack(originalBox);
202 }
203
204 SparseData::SparseData(const unsigned int dim, double min[], double max[])
205 {
206 // Creates a SparseData convering the range defined by min
207 // and max. For this it will create an empty Box for that
208 // range.
209 vector<double> minv(min,min+dim);
210 vector<double> maxv(max,max+dim);
211 Box originalBox(minv, maxv);
212 fList = new ProxyListBox();
213 fList->PushBack(originalBox);
214 }
215
217 { delete fList; }
218
219 unsigned int SparseData::NPoints() const
220 {
221 // Returns the number of points stored, including the 0 ones.
222 return fList->GetList().size();
223 }
224
225 unsigned int SparseData::NDim() const
226 {
227 // Returns the number of dimension of the SparseData object.
228 return fList->Begin()->GetMin().size();
229 }
230
231 void SparseData::Add(std::vector<double>& min, std::vector<double>& max,
232 const double content, const double error)
233 {
234 // Add a box to the stored ones. For that, it will look for
235 // the box that contains the new data and either replace it
236 // or updated it.
237
238 // Little box is the new Bin to be added
239 Box littleBox(min, max);
240 list<Box>::iterator it;
241 // So we look for the Bin already in the list that contains
242 // littleBox
243 it = std::find_if(fList->Begin(), fList->End(), BoxContainer(littleBox));
244 if ( it != fList->End() )
245// cout << "Found: " << *it << endl;
246 ;
247 else {
248 cout << "SparseData::Add -> FAILED! box not found! " << endl;
249 cout << littleBox << endl;
250 return; // Does not add the box, as it is part of the
251 // underflow/overflow bin
252 }
253 // If it happens to have a value, then we add the value,
254 if ( it->GetVal() )
255 it->AddVal( content );
256 else
257 {
258 // otherwise, we divide the container!
259 DivideBox(it->GetMin(), it->GetMax(),
260 littleBox.GetMin(), littleBox.GetMax(),
261 it->GetMin().size(), it->GetMin().size() - 1,
262 fList->GetList(), content, error );
263 // and remove it from the list
264 fList->Remove(it);
265 }
266 }
267
268 void SparseData::GetPoint(const unsigned int i,
269 std::vector<double>& min, std::vector<double>&max,
270 double& content, double& error)
271 {
272 // Get the point number i. This is a method to explore the
273 // data stored in the class.
274
275 unsigned int counter = 0;
276 list<Box>::iterator it = fList->Begin();
277 while ( it != fList->End() && counter != i ) {
278 ++it;
279 ++counter;
280 }
281
282 if ( (it == fList->End()) || (counter != i) )
283 throw std::out_of_range("SparseData::GetPoint");
284
285 min = it->GetMin();
286 max = it->GetMax();
287 content = it->GetVal();
288 error = it->GetError();
289 }
290
292 {
293 // Debug method to print a list with all the data stored.
294 copy(fList->Begin(), fList->End(), ostream_iterator<Box>(cout, "\n------\n"));
295 }
296
297
299 {
300 // Created the corresponding BinData
301
302 list<Box>::iterator it = fList->Begin();
303 const unsigned int dim = it->GetMin().size();
304
305 bd.Initialize(fList->GetList().size(), dim);
306 // Visit all the stored Boxes
307 for ( ; it != fList->End(); ++it )
308 {
309 vector<double> mid(dim);
310 // fill up the vector with the mid point of the Bin
311 for ( unsigned int i = 0; i < dim; ++i)
312 {
313 mid[i] = ((it->GetMax()[i] - it->GetMin()[i]) /2) + it->GetMin()[i];
314 }
315 // And store it into the BinData structure
316 bd.Add(&mid[0], it->GetVal(), it->GetError());
317 }
318 }
319
321 {
322 // Created the corresponding BinData as with the Integral
323 // option.
324
325 list<Box>::iterator it = fList->Begin();
326
327 bd.Initialize(fList->GetList().size(), it->GetMin().size());
328 // Visit all the stored Boxes
329 for ( ; it != fList->End(); ++it )
330 {
331 //Store the minimum value
332 bd.Add(&(it->GetMin()[0]), it->GetVal(), it->GetError());
333 //and the maximum
334 bd.AddBinUpEdge(&(it->GetMax()[0]));
335 }
336 }
337
339 {
340 // Created the corresponding BinData, but it does not include
341 // all the data with value equal to 0.
342
343 list<Box>::iterator it = fList->Begin();
344 const unsigned int dim = it->GetMin().size();
345
346 bd.Initialize(fList->GetList().size(), dim);
347 // Visit all the stored Boxes
348 for ( ; it != fList->End(); ++it )
349 {
350 // if the value is zero, jump to the next
351 if ( it->GetVal() == 0 ) continue;
352 vector<double> mid(dim);
353 // fill up the vector with the mid point of the Bin
354 for ( unsigned int i = 0; i < dim; ++i)
355 {
356 mid[i] = ((it->GetMax()[i] - it->GetMin()[i]) /2) + it->GetMin()[i];
357 }
358 // And store it into the BinData structure
359 bd.Add(&mid[0], it->GetVal(), it->GetError());
360 }
361 }
362
363 // Just for debugging pourposes
364 ostream& operator <<(ostream& os, const ROOT::Fit::Box& b)
365 {
366 os << "min: ";
367 copy(b.GetMin().begin(), b.GetMin().end(), ostream_iterator<double>(os, " "));
368 os << "max: ";
369 copy(b.GetMax().begin(), b.GetMax().end(), ostream_iterator<double>(os, " "));
370 os << "val: " << b.GetVal();
371
372 return os;
373 }
374 } // end namespace Fit
375
376} // end namespace ROOT
double
#define b(i)
Definition RSha256.hxx:100
size_t size(const MatrixT &matrix)
retrieve the size of a square matrix
void operator()(double value)
AreaComparer(vector< double >::iterator iter)
vector< double >::iterator fIter
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
bool operator()(const Box &b1)
BoxContainer(const Box &b)
vector< double > fMin
vector< double > fMax
const vector< double > & GetMax() const
void AddVal(const double value)
friend ostream & operator<<(ostream &os, const Box &b)
double GetVal() const
bool operator==(const Box &b)
Box(const vector< double > &min, const vector< double > &max, const double value=0.0, const double error=1.0)
const vector< double > & GetMin() const
double GetError() const
list< Box >::iterator Begin()
void PushBack(Box &box)
void Remove(list< Box >::iterator it)
list< Box > & GetList()
list< Box >::iterator End()
ProxyListBox * fList
Definition SparseData.h:61
void Add(std::vector< double > &min, std::vector< double > &max, const double content, const double error=1.0)
unsigned int NDim() const
void GetPoint(const unsigned int i, std::vector< double > &min, std::vector< double > &max, double &content, double &error)
unsigned int NPoints() const
void GetBinDataNoZeros(BinData &) const
void GetBinData(BinData &) const
SparseData(std::vector< double > &min, std::vector< double > &max)
void GetBinDataIntegral(BinData &) const
void box(Int_t pat, Double_t x1, Double_t y1, Double_t x2, Double_t y2)
Definition fillpatterns.C:1
const Int_t n
Definition legend1.C:16
TFitResultPtr Fit(FitObject *h1, TF1 *f1, Foption_t &option, const ROOT::Math::MinimizerOptions &moption, const char *goption, ROOT::Fit::DataRange &range)
Definition HFitImpl.cxx:133
void DivideBox(const vector< double > &min, const vector< double > &max, const vector< double > &bmin, const vector< double > &bmax, const unsigned int size, const unsigned int n, list< Box > &l, const double val, const double error)
tbb::task_arena is an alias of tbb::interface7::task_arena, which doesn't allow to forward declare tb...
auto * l
Definition textangle.C:4
REAL epsilon
Definition triangle.c:618