Logo ROOT   6.13/01
Reference Guide
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 
29 using namespace std;
30 
31 namespace 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
78  class BoxContainer
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.
119  class AreaComparer
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:
140  bool fThereIsArea;
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 
181  class ProxyListBox
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 
216  SparseData::~SparseData()
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 
291  void SparseData::PrintList() const
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 
298  void SparseData::GetBinData(BinData& bd) const
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 
320  void SparseData::GetBinDataIntegral(BinData& bd) const
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 
338  void SparseData::GetBinDataNoZeros(BinData& bd) const
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
you should not use this method at all Int_t Int_t Double_t Double_t Double_t Int_t mid
Definition: TRolke.cxx:630
Namespace for new ROOT classes and functions.
Definition: TFoamSampler.h:19
TMatrixT< Element > & Add(TMatrixT< Element > &target, Element scalar, const TMatrixT< Element > &source)
Modify addition: target += scalar * source.
Definition: TMatrixT.cxx:2933
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)
Definition: SparseData.cxx:152
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:627
std::ostream & operator<<(std::ostream &out, const TComplex &c)
Definition: TComplex.cxx:41
void Initialize(unsigned int newPoints, unsigned int dim=1, ErrorType err=kValueError)
Definition: BinData.cxx:349
Bool_t operator==(const TMatrixTBase< Element > &m1, const TMatrixTBase< Element > &m2)
Check to see if two matrices are identical.
VecExpr< UnaryOp< Fabs< T >, VecExpr< A, T, D >, T >, T, D > fabs(const VecExpr< A, T, D > &rhs)
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
void Add(double x, double y)
add one dim data with only coordinate and values
Definition: BinData.cxx:422
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