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  {
47  addBins(binedges);
48  }
49 
50 
56  Axis1D(const std::vector<BIN1D>& bins)
57  : _locked(false)
58  {
59  addBins(bins);
60  }
61 
62 
65  Axis1D(size_t nbins, double lower, double upper)
66  : _locked(false)
67  {
68  addBins(linspace(nbins, lower, upper));
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  {
80  addBins(bins);
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) {
359  _bins[i] += toAdd.bins().at(i);
360  }
361 
362  _dbn += toAdd._dbn;
363  _underflow += toAdd._underflow;
364  _overflow += toAdd._overflow;
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)
Add two axes together.
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
void addBin(double low, double high)
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
void addBins(const Bins &bins)
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