yoda is hosted by Hepforge, IPPP Durham
YODA - Yet more Objects for Data Analysis  1.7.2
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 
117  std::vector<double> xEdges() const {
118  std::vector<double> rtn(_binsearcher.edges().begin()+1, _binsearcher.edges().end()-1);
119  return rtn;
120  }
121 
123  BIN1D& bin(size_t index) {
124  if (index >= numBins()) throw RangeError("YODA::Histo1D: index out of range!");
125  return _bins[index];
126  }
127 
129  const BIN1D& bin(size_t index) const {
130  if (index >= numBins()) throw RangeError("YODA::Histo1D: index out of range!");
131  return _bins[index];
132  }
133 
135  ssize_t binIndexAt(double coord) const {
136  // Yes, this is robust even with an empty axis: there's always at least one outflow
137  return _indexes[_binsearcher.index(coord)];
138  }
139 
141  BIN1D& binAt(double x) {
142  const ssize_t index = binIndexAt(x);
143  if (index == -1) throw RangeError("There is no bin at the specified x");
144  return bin(index);
145  }
146 
148  const BIN1D& binAt(double x) const {
149  const ssize_t index = binIndexAt(x);
150  if (index == -1) throw RangeError("There is no bin at the specified x");
151  return bin(index);
152  }
153 
154 
156  const DBN& totalDbn() const {
157  return _dbn;
158  }
160  DBN& totalDbn() {
161  return _dbn;
162  }
164  void setTotalDbn(const DBN& dbn) {
165  _dbn = dbn;
166  }
167 
169  const DBN& underflow() const {
170  return _underflow;
171  }
173  DBN& underflow() {
174  return _underflow;
175  }
177  void setUnderflow(const DBN& dbn) {
178  _underflow = dbn;
179  }
180 
182  const DBN& overflow() const {
183  return _overflow;
184  }
186  DBN& overflow() {
187  return _overflow;
188  }
190  void setOverflow(const DBN& dbn) {
191  _overflow = dbn;
192  }
193 
195 
196 
198 
199 
201  void reset() {
202  _dbn.reset();
203  _underflow.reset();
204  _overflow.reset();
205  for (Bin& bin : _bins) bin.reset();
206  _locked = false;
207  }
208 
209 
211  void _setLock(bool locked) {
212  _locked = locked;
213  }
214 
215 
219  std::pair< std::vector<double>, std::vector<long> > _mk_edges_indexes(Bins& bins) const {
220  std::vector<double> edges; edges.reserve(bins.size()+1); // Nbins+1 edges
221  std::vector<long> indexes; edges.reserve(bins.size()+2); // Nbins + 2*outflows
222 
223  // Sort the bins
224  std::sort(bins.begin(), bins.end());
225 
226  // Keep a note of the last high edge
227  double last_high = -std::numeric_limits<double>::infinity();
228 
229  // Check for overlaps
230  for (size_t i = 0; i < bins.size(); ++i) {
231  Bin& currentBin = bins[i];
232  const double new_low = currentBin.xMin();
233  const double reldiff = (new_low - last_high) / currentBin.xWidth();
234  if (reldiff < -1e-3) { //< @note If there is a "large" negative gap (i.e. overlap), throw an exception
235  std::stringstream ss;
236  ss << "Bin edges overlap: " << last_high << " -> " << new_low;
237  throw RangeError(ss.str());
238  } else if (reldiff > 1e-3) { //< @note If there is a "large" positive gap, create a bin gap
239  indexes.push_back(-1); // First index will be for underflow
240  edges.push_back(new_low); // Will include first edge
241  }
242 
243  // Bins check that they are not zero or negative width. It's okay for
244  // them to throw an exception here, as we haven't changed anything yet.
245  indexes.push_back(i);
246  edges.push_back(currentBin.xMax());
247 
248  last_high = currentBin.xMax();
249  }
250  indexes.push_back(-1); // Overflow
251 
252  return std::make_pair(edges, indexes);
253  }
254 
255 
258  void mergeBins(size_t from, size_t to) {
259  // Correctness checking
260  if (from >= numBins())
261  throw RangeError("Initial merge index is out of range");
262  if (to >= numBins())
263  throw RangeError("Final merge index is out of range");
264  if (from > to)
265  throw RangeError("Final bin must be greater than or equal to initial bin");
266  if (_gapInRange(from, to))
267  throw RangeError("Bin ranges containing gaps cannot be merged");
268  if (from == to)
269  return; // nothing to be done
270 
271  Bin& b = bin(from);
272  for (size_t i = from + 1; i <= to; ++i)
273  b.merge(_bins[i]);
274  eraseBins(from+1, to);
275  }
276 
277 
283  void rebinBy(unsigned int n, size_t begin=0, size_t end=UINT_MAX) {
284  if (n < 1) throw UserError("Rebinning requested in groups of 0!");
285  for (size_t m = begin; m < end; ++m) {
286  if (m > numBins()) break;
287  const size_t myend = (m+n-1 < numBins()) ? m+n-1 : numBins()-1;
288  if (myend > m) {
289  mergeBins(m, myend);
290  end -= myend-m; //< reduce upper index by the number of removed bins
291  }
292  }
293  }
294 
296  void rebin(unsigned int n, size_t begin=0, size_t end=UINT_MAX) { rebinBy(n, begin, end); }
297 
298 
300  void rebinTo(const std::vector<double>& newedges) {
301  if (newedges.size() < 2)
302  throw UserError("Requested rebinning to an edge list which defines no bins");
303  const Utils::BinSearcher newbs(newedges);
304  const std::vector<double> eshared = newbs.shared_edges(_binsearcher);
305  if (eshared.size() != newbs.size())
306  throw BinningError("Requested rebinning to incompatible edges");
307  // std::cout << "Before merging" << std::endl;
308  // for (double x : _binsearcher.edges()) std::cout << x << std::endl;
309  // If the new min finite edge isn't the same, merge it into the underflow
310  // NB. Edge search match the *next* bin, so step back one unit... and note these are BinSearcher indices, i.e. i+1
311  if (!fuzzyEquals(xMin(), newedges.front())) {
312  const size_t kmatch = _binsearcher.index(newedges.front()) - 1;
313  mergeBins(0, kmatch-1);
314  _underflow += bin(0).dbn();
315  eraseBin(0);
316  }
317  // std::cout << "Merged start bins" << std::endl;
318  // for (double x : _binsearcher.edges()) std::cout << x << std::endl;
319  // Now the same for the overflow
320  if (!fuzzyEquals(xMax(), newedges.back())) {
321  const size_t kmatch = _binsearcher.index(newedges.back()) - 1;
322  // std::cout << newedges.back() << " -> " << kmatch << " .. " << _bins.size()-1 << " / " << numBins() << std::endl;
323  mergeBins(kmatch, _bins.size()-1);
324  _overflow += bin(_bins.size()-1).dbn();
325  eraseBin(_bins.size()-1);
326  }
327  // std::cout << "Merged end bins" << std::endl;
328  // for (double x : _binsearcher.edges()) std::cout << x << std::endl;
329  // Now merge the in-range bins
330  size_t jcurr = 0;
331  for (size_t i = 1; i < newedges.size(); ++i) { //< we already know that i=0 matches (until we support merging into overflows)
332  const size_t kmatch = _binsearcher.index(newedges.at(i)) - 1; //< Will match the *next* bin, so step back one unit... and note these are BinSearcher indices
333  assert(kmatch >= jcurr+1);
334  mergeBins(jcurr, kmatch-1);
335  jcurr += 1; //< The next bin to be merged, in the new numbering
336  }
337  // std::cout << "After merging" << std::endl;
338  // for (double x : _binsearcher.edges()) std::cout << x << std::endl;
339  }
340 
342  void rebin(const std::vector<double>& newedges) { rebinTo(newedges); }
343 
344 
346  void addBin(const Bin& b) {
348  Bins newBins(_bins);
349  newBins.push_back(b);
350  _updateAxis(newBins);
351  }
352 
353 
355  void addBin(double low, double high) {
356  addBin(Bin(low, high));
357  }
358 
359 
361  void addBins(const std::vector<double>& binedges) {
362  Bins newBins(_bins);
363  if (binedges.size() == 0) return;
364 
365  double low = binedges.front();
366  for (size_t i = 1; i < binedges.size(); ++i) {
367  const double high = binedges[i];
368  assert(high>low); // Make sure binedges are meaningful
369  newBins.push_back(Bin(low, high));
370  low = high;
371  }
372 
373  _updateAxis(newBins);
374  }
375 
376 
378  void addBins(const std::vector<std::pair<double, double> >& binpairs) {
379  // Make a copy of the current binning
380  Bins newBins(_bins);
381 
382  // Iterate over given bins
383  for (size_t i = 0; i < binpairs.size(); ++i) {
384  std::pair<double, double> b = binpairs[i];
385  newBins.push_back(Bin(b.first, b.second));
386  }
387  _updateAxis(newBins);
388  }
389 
390 
392  void addBins(const Bins& bins) {
393  Bins newBins(_bins);
394  for (const Bin& b : bins) newBins.push_back(b);
395  _updateAxis(newBins);
396  }
397 
398 
400  void eraseBin(const size_t i) {
401  // Might as well erase from the internal bins, as we can guarantee
402  // consistency.
403  if (i >= numBins())
404  throw RangeError("Bin index is out of range");
405 
406  const bool wasLocked = _locked;
407  _locked = false;
408  _bins.erase(_bins.begin() + i);
409  _updateAxis(_bins);
410  _locked = wasLocked;
411  }
412 
413 
415  void eraseBins(const size_t from, const size_t to) {
416  if (from >= numBins())
417  throw RangeError("Initial index out of range");
418  if (to >= numBins())
419  throw RangeError("Final index out of range");
420  if (from > to)
421  throw RangeError("Final index is less than initial index");
422 
423  const bool wasLocked = _locked;
424  _locked = false;
425  _bins.erase(_bins.begin() + from, _bins.begin() + to + 1);
426  _updateAxis(_bins);
427  _locked = wasLocked;
428  }
429 
430 
432  // @todo What if somebody passes in negative scalefactor? (test idea)
433  void scaleX(double scalefactor) {
434  _dbn.scaleX(scalefactor);
435  _underflow.scaleX(scalefactor);
436  _overflow.scaleX(scalefactor);
437  for (size_t i = 0; i < _bins.size(); ++i)
438  _bins[i].scaleX(scalefactor);
439  _updateAxis(_bins);
440  }
441 
442 
444  void scaleW(double scalefactor) {
445  _dbn.scaleW(scalefactor);
446  _underflow.scaleW(scalefactor);
447  _overflow.scaleW(scalefactor);
448  for (size_t i = 0; i < _bins.size(); ++i) _bins[i].scaleW(scalefactor);
449  }
450 
452 
453 
455 
456 
457  bool sameBinning(const Axis1D& other) const {
458  if (numBins() != other.numBins()) return false;
459  if (_indexes != other._indexes) return false;
460  return _binsearcher.same_edges(other._binsearcher);
461  }
462 
463  bool subsetBinning(const Axis1D& other) const {
464  const int ndiff = numBins() - other.numBins();
465  if (ndiff == 0) return sameBinning(other);
467  return !_binsearcher.shared_edges(other._binsearcher).empty();
468  }
469 
471 
472 
474 
475 
477  bool operator == (const Axis1D& other) const {
478  return sameBinning(other);
479  }
480 
481 
483  bool operator != (const Axis1D& other) const {
484  return ! operator == (other);
485  }
486 
487 
490  if (*this != toAdd) throw LogicError("YODA::Histo1D: Cannot add axes with different binnings.");
491 
492  for (size_t i = 0; i < _bins.size(); ++i) {
493  _bins[i] += toAdd.bins().at(i);
494  }
495 
496  _dbn += toAdd._dbn;
497  _underflow += toAdd._underflow;
498  _overflow += toAdd._overflow;
499  return *this;
500  }
501 
502 
505  if (*this != toSubtract) throw LogicError("YODA::Histo1D: Cannot add axes with different binnings.");
506 
507  for (size_t i = 0; i < _bins.size(); ++i) {
508  _bins[i] -= toSubtract.bins().at(i);
509  }
510 
511  _dbn -= toSubtract._dbn;
512  _underflow -= toSubtract._underflow;
513  _overflow -= toSubtract._overflow;
514  return *this;
515  }
516 
518 
519 
520  private:
521 
523  //
526  void _updateAxis(Bins& bins) {
527  // Ensure that axis is not locked
528  if (_locked) {
529  throw LockError("Attempting to update a locked 1D axis");
530  }
531 
532  // Get the new cuts and indexes (throws if overlaps), and set them on the searcher
533  const std::pair< std::vector<double>, std::vector<long> > es_is = _mk_edges_indexes(bins);
534  _binsearcher = Utils::BinSearcher(es_is.first);
535  _indexes = es_is.second;
536  _bins = bins;
537  }
538 
539 
542  bool _gapInRange(size_t ifrom, size_t ito) const {
543  if (ifrom == ito) return false;
544  assert(ito < numBins() && ifrom < ito);
545 
547  const size_t from_ix = _binsearcher.index(bin(ifrom).xMid());
548  const size_t to_ix = _binsearcher.index(bin(ito).xMid());
549  // std::cout << ifrom << " vs. " << from_ix << std::endl;
550  // std::cout << ito << " vs. " << to_ix << std::endl;
551 
552  for (size_t i = from_ix; i <= to_ix; i++)
553  // for (size_t i = ifrom; i <= ito; i++)
554  if (_indexes[i] == -1) return true;
555  return false;
556  }
557 
558 
559  private:
560 
562 
563 
565  Bins _bins;
566 
568  DBN _dbn;
569 
571  DBN _underflow;
572  DBN _overflow;
573 
574  // Binsearcher, for searching bins
575  Utils::BinSearcher _binsearcher;
576  // Mapping from binsearcher indices to bin indices (allowing gaps)
577  std::vector<long> _indexes;
578 
580  bool _locked;
581 
583 
584  };
585 
586 
588  template <typename BIN1D, typename DBN>
589  inline Axis1D<BIN1D,DBN> operator + (const Axis1D<BIN1D,DBN>& first, const Axis1D<BIN1D,DBN>& second) {
590  Axis1D<BIN1D,DBN> tmp = first;
591  tmp += second;
592  return tmp;
593  }
594 
596  template <typename BIN1D, typename DBN>
597  inline Axis1D<BIN1D,DBN> operator - (const Axis1D<BIN1D,DBN>& first, const Axis1D<BIN1D,DBN>& second) {
598  Axis1D<BIN1D,DBN> tmp = first;
599  tmp -= second;
600  return tmp;
601  }
602 
603 
604 }
605 
606 #endif
BIN1D Bin
Typedefs.
Definition: Axis1D.h:27
double xMin() const
Return the lowest-value bin edge on the axis.
Definition: Axis1D.h:102
Axis1D< BIN1D, DBN > & operator+=(const Axis1D< BIN1D, DBN > &toAdd)
Add two axes together.
Definition: Axis1D.h:489
Axis1D(size_t nbins, double lower, double upper)
Definition: Axis1D.h:65
void addBin(double low, double high)
Add a bin, providing its low and high edge.
Definition: Axis1D.h:355
bool operator!=(const Axis1D &other) const
Check if the binning of two of the Axis is different.
Definition: Axis1D.h:483
void eraseBin(const size_t i)
Remove a bin.
Definition: Axis1D.h:400
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:141
DBN & underflow()
Return underflow (non-const)
Definition: Axis1D.h:173
Axis1D< BIN1D, DBN > operator-(const Axis1D< BIN1D, DBN > &first, const Axis1D< BIN1D, DBN > &second)
Subtract the statistics on two axis.
Definition: Axis1D.h:597
void addBins(const std::vector< std::pair< double, double > > &binpairs)
Add a list of bins as pairs of lowEdge, highEdge.
Definition: Axis1D.h:378
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:72
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:252
void eraseBins(const size_t from, const size_t to)
Remove a bin range.
Definition: Axis1D.h:415
DBN & overflow()
Return overflow (non-const)
Definition: Axis1D.h:186
BIN1D & bin(size_t index)
Return a bin at a given index (non-const)
Definition: Axis1D.h:123
const BIN1D & binAt(double x) const
Return a bin at a given coordinate (const)
Definition: Axis1D.h:148
void scaleX(double scalefactor)
Scale the size of an axis by a factor.
Definition: Axis1D.h:433
const DBN & underflow() const
Return underflow (const)
Definition: Axis1D.h:169
Error for problems introduced outside YODA, to put it nicely.
Definition: Exceptions.h:100
void rebinBy(unsigned int n, size_t begin=0, size_t end=UINT_MAX)
Merge every group of n bins, from start to end inclusive.
Definition: Axis1D.h:283
void addBins(const std::vector< double > &binedges)
Add a contiguous set of bins to an axis, via their list of edges.
Definition: Axis1D.h:361
const DBN & overflow() const
Return overflow (const)
Definition: Axis1D.h:182
Axis1D< BIN1D, DBN > & operator-=(const Axis1D< BIN1D, DBN > &toSubtract)
Subtract two axes.
Definition: Axis1D.h:504
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
Error for general binning problems.
Definition: Exceptions.h:27
Axis1D()
Empty constructor.
Definition: Axis1D.h:38
const DBN & totalDbn() const
Return the total distribution (const)
Definition: Axis1D.h:156
void rebin(const std::vector< double > &newedges)
Overloaded alias for rebinTo.
Definition: Axis1D.h:342
Axis1D< BIN1D, DBN > operator+(const Axis1D< BIN1D, DBN > &first, const Axis1D< BIN1D, DBN > &second)
Add the statistics on two axes.
Definition: Axis1D.h:589
Bins & bins()
Return a vector of bins (non-const)
Definition: Axis1D.h:97
void mergeBins(size_t from, size_t to)
Definition: Axis1D.h:258
std::vector< double > xEdges() const
Definition: Axis1D.h:117
size_t numBins() const
Get the number of bins on the axis.
Definition: Axis1D.h:87
void reset()
Reset all the bin statistics on the axis.
Definition: Axis1D.h:201
void addBins(const Bins &bins)
Add a list of Bin objects.
Definition: Axis1D.h:392
bool sameBinning(const Axis1D &other) const
Definition: Axis1D.h:457
const Bins & bins() const
Return a vector of bins (const)
Definition: Axis1D.h:92
std::vector< Bin > Bins
A vector containing 1D bins. Not used for searching.
Definition: Axis1D.h:30
void addBin(const Bin &b)
Add a bin, passed explicitly.
Definition: Axis1D.h:346
void setUnderflow(const DBN &dbn)
Set the underflow distribution: CAREFUL!
Definition: Axis1D.h:177
Axis1D(const std::vector< BIN1D > &bins)
Constructor accepting a vector of bins.
Definition: Axis1D.h:56
double xMax() const
Return the highest-value bin edge on the axis.
Definition: Axis1D.h:108
void rebinTo(const std::vector< double > &newedges)
Rebin to the given list of bin edges.
Definition: Axis1D.h:300
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
void setTotalDbn(const DBN &dbn)
Set the total distribution: CAREFUL!
Definition: Axis1D.h:164
const BIN1D & bin(size_t index) const
Return a bin at a given index (const)
Definition: Axis1D.h:129
DBN & totalDbn()
Return the total distribution (non-const)
Definition: Axis1D.h:160
void rebin(unsigned int n, size_t begin=0, size_t end=UINT_MAX)
Overloaded alias for rebinBy.
Definition: Axis1D.h:296
bool subsetBinning(const Axis1D &other) const
Definition: Axis1D.h:463
ssize_t binIndexAt(double coord) const
Returns an index of a bin at a given coord, -1 if no bin matches.
Definition: Axis1D.h:135
void setOverflow(const DBN &dbn)
Set the overflow distribution: CAREFUL!
Definition: Axis1D.h:190
bool operator==(const Axis1D &other) const
Check if two of the Axis have the same binning, within numeric tolerance.
Definition: Axis1D.h:477
void scaleW(double scalefactor)
Scale the amount of fills by a factor.
Definition: Axis1D.h:444
1D bin container
Definition: Axis1D.h:20