YODA is hosted by Hepforge, IPPP Durham
 YODA - Yet more Objects for Data Analysis  1.3.1
Axis1D.h
Go to the documentation of this file.
1 #ifndef YODA_Axis1D_h
2 #define YODA_Axis1D_h
3
4 #include "YODA/AnalysisObject.h"
5 #include "YODA/Exceptions.h"
6 #include "YODA/Bin.h"
7 #include "YODA/Utils/MathUtils.h"
9 #include <limits>
10 #include <string>
11
12 namespace YODA {
13
14
19  template <typename BIN1D, typename DBN>
20  class Axis1D {
21  public:
22
24
25
27  typedef BIN1D Bin;
28
30  typedef typename std::vector<Bin> Bins;
31
33
35
36
39  : _locked(false)
40  { }
41
42
44  Axis1D(const std::vector<double>& binedges)
45  : _locked(false)
46  {
48  }
49
50
56  Axis1D(const std::vector<BIN1D>& bins)
57  : _locked(false)
58  {
60  }
61
62
65  Axis1D(size_t nbins, double lower, double upper)
66  : _locked(false)
67  {
69  }
70
71
77  Axis1D(const Bins& bins, const DBN& dbn_tot, const DBN& dbn_uflow, const DBN& dbn_oflow)
78  : _dbn(dbn_tot), _underflow(dbn_uflow), _overflow(dbn_oflow), _locked(false)
79  {
81  }
82
84
85
87  size_t numBins() const {
88  return bins().size();
89  }
90
92  const Bins& bins() const {
93  return _bins;
94  }
95
97  Bins& bins() {
98  return _bins;
99  }
100
102  double xMin() const {
103  if (numBins() == 0) throw RangeError("This axis contains no bins and so has no defined range");
104  return bins().front().xMin();
105  }
106
108  double xMax() const {
109  if (numBins() == 0) throw RangeError("This axis contains no bins and so has no defined range");
110  return bins().back().xMax();
111  }
112
114  BIN1D& bin(size_t index) {
115  if (index >= numBins()) throw RangeError("YODA::Histo1D: index out of range!");
116  return _bins[index];
117  }
118
120  const BIN1D& bin(size_t index) const {
121  if (index >= numBins()) throw RangeError("YODA::Histo1D: index out of range!");
122  return _bins[index];
123  }
124
126  ssize_t binIndexAt(double coord) const {
127  // Yes, this is robust even with an empty axis: there's always at least one outflow
128  return _indexes[_binsearcher.index(coord)];
129  }
130
132  BIN1D& binAt(double x) {
133  const ssize_t index = binIndexAt(x);
134  if (index == -1) throw RangeError("There is no bin at the specified x");
135  return bin(index);
136  }
137
139  const BIN1D& binAt(double x) const {
140  const ssize_t index = binIndexAt(x);
141  if (index == -1) throw RangeError("There is no bin at the specified x");
142  return bin(index);
143  }
144
145
147  const DBN& totalDbn() const {
148  return _dbn;
149  }
150
152  DBN& totalDbn() {
153  return _dbn;
154  }
155
157  const DBN& underflow() const {
158  return _underflow;
159  }
160
162  DBN& underflow() {
163  return _underflow;
164  }
165
167  const DBN& overflow() const {
168  return _overflow;
169  }
170
172  DBN& overflow() {
173  return _overflow;
174  }
175
176
178
179
181
182
184  void reset() {
185  _dbn.reset();
186  _underflow.reset();
187  _overflow.reset();
188  BOOST_FOREACH(Bin& bin, _bins) bin.reset();
189  _locked = false;
190  }
191
192
194  void _setLock(bool locked) {
195  _locked = locked;
196  }
197
198
201  void mergeBins(size_t from, size_t to) {
202  // Correctness checking
203  if (from >= numBins())
204  throw RangeError("Initial merge index is out of range");
205  if (to >= numBins())
206  throw RangeError("Final merge index is out of range");
207  if (from > to)
208  throw RangeError("Final bin must be greater than initial bin");
209  if (_gapInRange(from, to))
210  throw RangeError("Bin ranges containing gaps cannot be merged");
211
212  Bin& b = bin(from);
213  for (size_t i = from + 1; i <= to; i++)
214  b.merge(_bins[i]);
215  eraseBins(from+1, to);
216  }
217
223  void rebin(size_t n) {
224  for (size_t m = 0; m < numBins(); m++) {
225  const size_t end = (m + n - 1 < numBins()) ? m + n -1 : numBins() - 1;
226  if (end > m) mergeBins(m, end);
227  }
228  }
229
231  void addBin(double low, double high) {
232  Bins newBins(_bins);
233  newBins.push_back(Bin(low, high));
234  _updateAxis(newBins);
235  }
236
237
239  void addBins(const std::vector<double>& binedges) {
240  Bins newBins(_bins);
241  if (binedges.size() == 0)
242  return;
243
244  double low = binedges.front();
245
246  for (size_t i = 1; i < binedges.size(); i++) {
247  double high = binedges[i];
248  newBins.push_back(Bin(low, high));
249  low = high;
250  }
251
252  _updateAxis(newBins);
253  }
254
255  // (The following proposed method would drastically simplify the
256  // Python interfaces handling of bin adding -- it would also
257  // simplify some of the hurdles of bins not having default
258  // constructors)
259
261  void addBins(const std::vector<std::pair<double, double> >& binpairs) {
262  // Make a copy of the current binning
263  Bins newBins(_bins);
264
265  // Iterate over given bins
266  for (size_t i = 0; i < binpairs.size(); i++) {
267  std::pair<double, double> b = binpairs[i];
268  newBins.push_back(Bin(b.first, b.second));
269  }
270  _updateAxis(newBins);
271  }
272
273
274  void addBins(const Bins& bins) {
275  Bins newBins(_bins);
276  BOOST_FOREACH(const Bin& b, bins) newBins.push_back(b);
277  _updateAxis(newBins);
278  }
279
281  void eraseBin(const size_t i) {
282  // Might as well erase from the internal bins, as we can guarantee
283  // consistency.
284  if (i >= numBins())
285  throw RangeError("Bin index is out of range");
286
287  const bool wasLocked = _locked;
288  _locked = false;
289  _bins.erase(_bins.begin() + i);
290  _updateAxis(_bins);
291  _locked = wasLocked;
292  }
293
295  void eraseBins(const size_t from, const size_t to) {
296  if (from >= numBins())
297  throw RangeError("Initial index out of range");
298  if (to >= numBins())
299  throw RangeError("Final index out of range");
300  if (from > to)
301  throw RangeError("Final index is less than initial index");
302
303  const bool wasLocked = _locked;
304  _locked = false;
305  _bins.erase(_bins.begin() + from, _bins.begin() + to + 1);
306  _updateAxis(_bins);
307  _locked = wasLocked;
308  }
309
311  // @todo What if somebody passes in negative scalefactor? (test idea)
312  void scaleX(double scalefactor) {
313  _dbn.scaleX(scalefactor);
314  _underflow.scaleX(scalefactor);
315  _overflow.scaleX(scalefactor);
316  for (size_t i = 0; i < _bins.size(); ++i)
317  _bins[i].scaleX(scalefactor);
318  _updateAxis(_bins);
319  }
320
322  void scaleW(double scalefactor) {
323  _dbn.scaleW(scalefactor);
324  _underflow.scaleW(scalefactor);
325  _overflow.scaleW(scalefactor);
326  for (size_t i = 0; i < _bins.size(); ++i) _bins[i].scaleW(scalefactor);
327  }
328
330
331
333
334
336  bool operator == (const Axis1D& other) const {
337  if (numBins() != other.numBins())
338  return false;
339  for (size_t i = 0; i < numBins(); i++)
340  if (!(fuzzyEquals(bin(i).xMin(), other.bin(i).xMin()) &&
341  fuzzyEquals(bin(i).xMax(), other.bin(i).xMax())))
342  return false;
343
344  return true;
345  }
346
347
349  bool operator != (const Axis1D& other) const {
350  return ! operator == (other);
351  }
352
353
356  if (*this != toAdd) throw LogicError("YODA::Histo1D: Cannot add axes with different binnings.");
357
358  for (size_t i = 0; i < _bins.size(); ++i) {
360  }
361
365  return *this;
366  }
367
368
371  if (*this != toSubtract) throw LogicError("YODA::Histo1D: Cannot add axes with different binnings.");
372
373  for (size_t i = 0; i < _bins.size(); ++i) {
374  _bins[i] -= toSubtract.bins().at(i);
375  }
376
377  _dbn -= toSubtract._dbn;
378  _underflow -= toSubtract._underflow;
379  _overflow -= toSubtract._overflow;
380  return *this;
381  }
382
384
385
386  private:
387
389  //
392  void _updateAxis(Bins& bins) {
393  // Ensure that axis is not locked
394  if (_locked) {
395  throw LockError("Attempting to update a locked axis");
396  }
397  // Define the new cuts and indexes
398  std::vector<double> edges; edges.reserve(bins.size()+1); // Nbins+1 edges
399  std::vector<long> indexes; edges.reserve(bins.size()+2); // Nbins + 2 outflows
400
401  // Sort the bins
402  std::sort(bins.begin(), bins.end());
403
404  // Keep a note of the last high edge
405  double last_high = -std::numeric_limits<double>::infinity();
406
407  // Check for overlaps
408  for (size_t i = 0; i < bins.size(); ++i) {
409  Bin& currentBin = bins[i];
410  const double new_low = currentBin.xMin();
411  const double reldiff = (new_low - last_high) / currentBin.xWidth();
412  if (reldiff < -1e-3) { //< @note If there is a "large" negative gap (i.e. overlap), throw an exception
413  std::stringstream ss;
414  ss << "Bin edges overlap: " << last_high << " -> " << new_low;
415  throw RangeError(ss.str());
416  } else if (reldiff > 1e-3) { //< @note If there is a "large" positive gap, create a bin gap
417  indexes.push_back(-1); // First index will be for underflow
418  edges.push_back(new_low); // Will include first edge
419  }
420
421  // Bins check that they are not zero or negative width. It's okay for
422  // them to throw an exception here, as we haven't changed anything yet.
423  indexes.push_back(i);
424  edges.push_back(currentBin.xMax());
425
426  last_high = currentBin.xMax();
427  }
428  indexes.push_back(-1); // Overflow
429
430  // Everything was okay, so let's make our changes
431  _binsearcher = Utils::BinSearcher(edges);
432  _indexes = indexes;
433  _bins = bins;
434  }
435
436
439  bool _gapInRange(size_t ifrom, size_t ito) const {
440  if (ifrom == ito) return false;
441  assert(ito < numBins() && ifrom < ito);
442
444  const size_t from_ix = _binsearcher.index(bin(ifrom).xMid());
445  const size_t to_ix = _binsearcher.index(bin(ito).xMid());
446  std::cout << ifrom << " vs. " << from_ix << std::endl;
447  std::cout << ito << " vs. " << to_ix << std::endl;
448
449  for (size_t i = from_ix; i <= to_ix; i++)
450  // for (size_t i = ifrom; i <= ito; i++)
451  if (_indexes[i] == -1) return true;
452  return false;
453  }
454
455
456  private:
457
459
460
462  Bins _bins;
463
465  DBN _dbn;
466
468  DBN _underflow;
469  DBN _overflow;
470
471  // Binsearcher, for searching bins
472  Utils::BinSearcher _binsearcher;
473  // Mapping from binsearcher indices to bin indices (allowing gaps)
474  std::vector<long> _indexes;
475
477  bool _locked;
478
480
481  };
482
483
485  template <typename BIN1D, typename DBN>
486  inline Axis1D<BIN1D,DBN> operator + (const Axis1D<BIN1D,DBN>& first, const Axis1D<BIN1D,DBN>& second) {
487  Axis1D<BIN1D,DBN> tmp = first;
488  tmp += second;
489  return tmp;
490  }
491
493  template <typename BIN1D, typename DBN>
494  inline Axis1D<BIN1D,DBN> operator - (const Axis1D<BIN1D,DBN>& first, const Axis1D<BIN1D,DBN>& second) {
495  Axis1D<BIN1D,DBN> tmp = first;
496  tmp -= second;
497  return tmp;
498  }
499
500 }
501
502 #endif
BIN1D Bin
Typedefs.
Definition: Axis1D.h:27
size_t numBins() const
Get the number of bins on the axis.
Definition: Axis1D.h:87
Axis1D< BIN1D, DBN > & operator+=(const Axis1D< BIN1D, DBN > &toAdd)
Definition: Axis1D.h:355
Axis1D(size_t nbins, double lower, double upper)
Definition: Axis1D.h:65
const BIN1D & binAt(double x) const
Return a bin at a given coordinate (const)
Definition: Axis1D.h:139
void rebin(size_t n)
Merge every group of n bins, starting from the LHS.
Definition: Axis1D.h:223
Add a bin, providing its low and high edge.
Definition: Axis1D.h:231
void eraseBin(const size_t i)
Remove a bin.
Definition: Axis1D.h:281
Axis1D(const std::vector< double > &binedges)
Constructor accepting a list of bin edges.
Definition: Axis1D.h:44
BIN1D & binAt(double x)
Return a bin at a given coordinate (non-const)
Definition: Axis1D.h:132
DBN & underflow()
Return underflow (non-const)
Definition: Axis1D.h:162
Axis1D< BIN1D, DBN > operator-(const Axis1D< BIN1D, DBN > &first, const Axis1D< BIN1D, DBN > &second)
Subtract the statistics on two axis.
Definition: Axis1D.h:494
const DBN & overflow() const
Return overflow (const)
Definition: Axis1D.h:167
void addBins(const std::vector< std::pair< double, double > > &binpairs)
Add a list of bins as pairs of lowEdge, highEdge.
Definition: Axis1D.h:261
Error for e.g. use of invalid bin ranges.
Definition: Exceptions.h:34
bool fuzzyEquals(double a, double b, double tolerance=1E-5)
Compare two floating point numbers for equality with a degree of fuzziness.
Definition: MathUtils.h:68
std::vector< double > linspace(size_t nbins, double start, double end, bool include_end=true)
Make a list of nbins + 1 values equally spaced between start and end inclusive.
Definition: MathUtils.h:248
void eraseBins(const size_t from, const size_t to)
Remove a bin range.
Definition: Axis1D.h:295
DBN & overflow()
Return overflow (non-const)
Definition: Axis1D.h:172
BIN1D & bin(size_t index)
Return a bin at a given index (non-const)
Definition: Axis1D.h:114
const BIN1D & bin(size_t index) const
Return a bin at a given index (const)
Definition: Axis1D.h:120
bool operator==(const Axis1D &other) const
Check if two of the Axis have the same binning.
Definition: Axis1D.h:336
double xMin() const
Return the lowest-value bin edge on the axis.
Definition: Axis1D.h:102
double xMax() const
Return the highest-value bin edge on the axis.
Definition: Axis1D.h:108
void scaleX(double scalefactor)
Scale the size of an axis by a factor.
Definition: Axis1D.h:312
const DBN & underflow() const
Return underflow (const)
Definition: Axis1D.h:157
void addBins(const std::vector< double > &binedges)
Add a contiguous set of bins to an axis, via their list of edges.
Definition: Axis1D.h:239
ssize_t binIndexAt(double coord) const
Returns an index of a bin at a given coord, -1 if no bin matches.
Definition: Axis1D.h:126
Axis1D< BIN1D, DBN > & operator-=(const Axis1D< BIN1D, DBN > &toSubtract)
Subtract two axes.
Definition: Axis1D.h:370
Error for modification of a data object where filling has already begun.
Definition: Exceptions.h:41
Error for places where it should not have been possible to get to!
Definition: Exceptions.h:55
const Bins & bins() const
Return a vector of bins (const)
Definition: Axis1D.h:92
Axis1D()
Empty constructor.
Definition: Axis1D.h:38
Axis1D< BIN1D, DBN > operator+(const Axis1D< BIN1D, DBN > &first, const Axis1D< BIN1D, DBN > &second)
Add the statistics on two axes.
Definition: Axis1D.h:486
Bins & bins()
Return a vector of bins (non-const)
Definition: Axis1D.h:97
void mergeBins(size_t from, size_t to)
Definition: Axis1D.h:201
void reset()
Reset all the bin statistics on the axis.
Definition: Axis1D.h:184
const DBN & totalDbn() const
Return the total distribution (const)
Definition: Axis1D.h:147
Definition: Axis1D.h:274
std::vector< Bin > Bins
A vector containing 1D bins. Not used for searching.
Definition: Axis1D.h:30
Axis1D(const std::vector< BIN1D > &bins)
Constructor accepting a vector of bins.
Definition: Axis1D.h:56
Axis1D(const Bins &bins, const DBN &dbn_tot, const DBN &dbn_uflow, const DBN &dbn_oflow)
Constructor accepting a vector of bins.
Definition: Axis1D.h:77
DBN & totalDbn()
Return the total distribution (non-const)
Definition: Axis1D.h:152
bool operator!=(const Axis1D &other) const
Check if the binning of two of the Axis is different.
Definition: Axis1D.h:349
void scaleW(double scalefactor)
Scale the amount of fills by a factor.
Definition: Axis1D.h:322
1D bin container
Definition: Axis1D.h:20