|
[yoda-svn] r329 - in trunk: . include/YODA srcblackhole at projects.hepforge.org blackhole at projects.hepforge.orgTue Aug 23 01:54:27 BST 2011
Author: buckley Date: Tue Aug 23 01:54:26 2011 New Revision: 329 Log: Adding first implementation of 1D bin merging to Axis1D. Untested. NB. Remember to add entries to the ChangeLog for significant commits. Modified: trunk/ChangeLog trunk/TODO trunk/include/YODA/Axis1D.h trunk/include/YODA/Bin1D.h trunk/include/YODA/Histo1D.h trunk/src/Histo1D.cc Modified: trunk/ChangeLog ============================================================================== --- trunk/ChangeLog Mon Aug 22 18:01:34 2011 (r328) +++ trunk/ChangeLog Tue Aug 23 01:54:26 2011 (r329) @@ -1,3 +1,7 @@ +2011-08-23 Andy Buckley <andy at insectnation.org> + + * Adding first implementation of 1D bin merging to Axis1D. + 2011-08-22 Andy Buckley <andy at insectnation.org> * Adding copy constructors and assignment operators to Histo1D, Modified: trunk/TODO ============================================================================== --- trunk/TODO Mon Aug 22 18:01:34 2011 (r328) +++ trunk/TODO Tue Aug 23 01:54:26 2011 (r329) @@ -4,9 +4,10 @@ NOW * Bin division: add quadratic / binomial error treatment enum option (AB) - AB: done... needs to be checked. Add a test that y +- 1sigma in [0, 1] + AB: done... needs to be checked. Add a test that y +- 1 sigma in [0, 1] * Rebinning in 1D: merges of n adjacent bins (AB) + AB: Done on Axis1D... add interface to Histo1D and Profile1D. * Add copy assignment to Point/Scatter3D and Dbn1D/2D. Modified: trunk/include/YODA/Axis1D.h ============================================================================== --- trunk/include/YODA/Axis1D.h Mon Aug 22 18:01:34 2011 (r328) +++ trunk/include/YODA/Axis1D.h Tue Aug 23 01:54:26 2011 (r329) @@ -9,58 +9,71 @@ #include <string> namespace YODA { + + /// @brief 1D bin container - /// This class handles most of the low-level operations on the bins - /// arranged in 1D line. + /// + /// This class handles most of the low-level operations on an axis of bins + /// arranged in a 1D line. template <typename BIN1D, typename DBN> class Axis1D { - public: + public: + + /// Convenience typedefs + typedef BIN1D Bin; - /// Convinience typedefs - typedef BIN1D Bin; + /// A vector containing 1D bins. It is not used for searching, + /// only for bins storage. + typedef typename std::vector<BIN1D> Bins; + + /// A standard edge container. First number defines location + /// of the edge, the second one a bin nuber the edge is a member of. + typedef typename std::pair<double, size_t> Edge; - /// A vector containing 1D bins. It is not used for searching, - /// only for bins storage. - typedef typename std::vector<BIN1D> Bins; - - /// A standard edge container. First number defines location - /// of the edge, the second one a bin nuber the edge is a member of. - typedef typename std::pair<double, size_t> Edge; /// @name Constructors //@{ + /// Empty constructor Axis1D() { } + /// Constructor accepting a list of bin edges Axis1D(const std::vector<double>& binedges) { _mkAxis(binedges); } + /// Constructor with histogram limits, number of bins and a bin distribution enum Axis1D(size_t nbins, double lower, double upper) { _mkAxis(linspace(lower, upper, nbins)); } + /// @brief Constructor accepting a vector of bins. + /// /// Note that not only dimensions of these bins will be set, - /// all the contents of the bins will be copied across, including + /// all the contents of the bins will be copied across, including /// the statistics Axis1D(const std::vector<BIN1D>& bins) { _mkAxis(bins); } + /// State-setting constructor Axis1D(const Bins& bins, const DBN& dbn_tot, const DBN& dbn_uflow, const DBN& dbn_oflow) : _dbn(dbn_tot), _underflow(dbn_uflow), _overflow(dbn_oflow) { _mkAxis(bins); } + //@} - + + /// @name Statistics accessor functions //@{ + /// Get the number of bins on the axis size_t numBins() const { return _bins.size(); @@ -70,7 +83,7 @@ Bins& bins() { return _bins; } - + /// Return a vector of bins (const) const Bins& bins() const { return _bins; @@ -95,13 +108,13 @@ if(index >= numBins()) throw RangeError("YODA::Histo1D: index out of range!"); return _bins[index]; } - + /// Return a bin at a given index (const) const BIN1D& bin (size_t index) const { if(index >= numBins()) throw RangeError("YODA::Histo1D: index out of range!"); return _bins[index]; } - + /// Return a bin at a given coordinate (non-const) BIN1D& binByCoord(double x) { int index = getBinIndex(x); @@ -115,7 +128,7 @@ if(index == -1) throw RangeError("There is no bin at the specified x"); return bin(index); } - + /// Return the total distribution (non-const) DBN& totalDbn() { return _dbn; @@ -145,102 +158,131 @@ const DBN& overflow() const { return _overflow; } + //@} + /// @name Modifiers and helpers //@{ + /// Returns an index of a bin at a given coord, -1 if no bin. int getBinIndex(double coord) const { - /// This is NOT a magic number, it is merely a 'trick' to yield a - /// correct answer in the case of coord pointing exactly on an edge. - /// This will never change the answer, since it is *much* smaller than - /// the fuzzyEquals() tolerance! - coord += 0.00000000001; - - /// Search for an edge + // This is NOT a magic number, it is merely a 'trick' to yield a + // correct answer in the case of coord pointing exactly on an edge. + // This will never change the answer, since it is *much* smaller than + // the fuzzyEquals() tolerance! + coord += 0.00000000001; + + // Search for an edge size_t index = _binaryS(coord, 0, _binHashSparse.size()); - /// If lower bound is a last edge, it means that we are 'above' the - /// axis and no bin canbe found in such case. Therefore, announce it. - if(index == _binHashSparse.size() - 1) return -1; - - /// If on the left to the point in consideration there is an edge that is - /// a member of the same bin as the one on the right, it means that our point - /// is inside a bin. In such case, announce it providing the index of the - /// bin in question. - if(_binHashSparse[index].second == _binHashSparse[index+1].second) - { + // If lower bound is a last edge, it means that we are 'above' the + // axis and no bin canbe found in such case. Therefore, announce it. + if (index == _binHashSparse.size() - 1) return -1; + + // If on the left to the point in consideration there is an edge that is + // a member of the same bin as the one on the right, it means that our point + // is inside a bin. In such case, announce it providing the index of the + // bin in question. + if (_binHashSparse[index].second == _binHashSparse[index+1].second) { return _binHashSparse[index].second; } - /// If we are inside an axis, but not inside a bin, it means that we must - /// be trying to fill an empty space. + // If we are inside an axis, but not inside a bin, it means that we must + // be trying to fill an empty space. return -1; } + /// Reset all the statistics on the Axis void reset() { _dbn.reset(); _underflow.reset(); _overflow.reset(); - for(size_t i = 0; i < _bins.size(); ++i) _bins[i].reset(); + for (size_t i = 0; i < _bins.size(); ++i) _bins[i].reset(); } - void rebin(int factor) { - throw std::runtime_error("Implement!"); + + /// Merge every group of n bins, starting from the LHS + void rebin(int n) { + int m = 0; + while (m < _bins.size()) { + size_t end = (m + n - 1 < _bins.size()) ? m + n -1 : _bins.size() - 1; + mergeBins(m, end); + m += 1; + } } - + + + /// Merge together the bin range with indices from @a from to @a to, inclusive void mergeBins(size_t from, size_t to) { - /// Correctness checking - if(from < 0 || from >= numBins()) throw ("First index is out of range!"); - if(to < 0 || to >= numBins()) throw ("Second index is out of range!"); - if(_bins[from].xMin() > _bins[to].xMin()) throw RangeError("The starting bin is greater than ending bin!"); - throw std::runtime_error("Implement!"); - //@todo Implement! + // Correctness checking + if (from < 0 || from >= numBins()) throw ("First index is out of range!"); + if (to < 0 || to >= numBins()) throw ("Second index is out of range!"); + if (_bins[from].xMin() > _bins[to].xMin()) throw RangeError("The starting bin is greater than ending bin!"); + BIN1D& b = _bins[from]; + for (size_t i = from+1; i <= to; ++i) { + b.merge(_bins[i]); + } + eraseBins(from+1, to); } + /// Add a bin, providing its low and high edge - void addBin(double& from, double& to) { + void addBin(double from, double to) { std::vector<double> binedges; binedges.push_back(from); binedges.push_back(to); - _mkAxis(binedges); } - /// Add a set of bins to an axis, providing their boundaries + + /// Add a set of bins to an axis, providing their boundaries. + /// /// Notice that no empty space will be created. - void addBin(std::vector<double>& binedges) { + void addBin(const std::vector<double>& binedges) { _mkAxis(binedges); } - /// Add a set of bins specifying start/end of each - void addBin(std::vector<std::pair<double,double> > edges) { - for(size_t i = 0; i < edges.size(); ++i){ + + /// Add a set of bins specifying start/end of each + void addBin(const std::vector<std::pair<double,double> > edges) { + for (size_t i = 0; i < edges.size(); ++i) { std::vector<double> temp; temp.push_back(edges[i].first); temp.push_back(edges[i].second); - _mkAxis(temp); } } + /// Remove a bin void eraseBin(size_t index) { - if(index >= _bins.size()) throw RangeError("Index out of range!"); + if (index >= _bins.size()) throw RangeError("Index out of range!"); _bins.erase(_bins.begin() + index); _mkBinHash(); } - + + + /// Remove a bin range + void eraseBins(size_t from, size_t to) { + if (from < 0 || from >= numBins()) throw ("First index is out of range!"); + if (to < 0 || to >= numBins()) throw ("Second index is out of range!"); + if (_bins[from].xMin() > _bins[to].xMin()) throw RangeError("The starting bin is greater than ending bin!"); + _bins.erase(_bins.begin() + from, _bins.begin() + to + 1); + _mkBinHash(); + } + + /// Scale the size of an axis by a factor void scaleX(double scalefactor) { _dbn.scaleX(scalefactor); _underflow.scaleX(scalefactor); _overflow.scaleW(scalefactor); - for(size_t i = 0; i < _bins.size(); ++i) _bins[i].scaleX(scalefactor); - + for (size_t i = 0; i < _bins.size(); ++i) _bins[i].scaleX(scalefactor); _mkBinHash(); } - + + /// Scale the amount of fills by a factor void scaleW(double scalefactor) { _dbn.scaleW(scalefactor); @@ -248,8 +290,10 @@ _overflow.scaleW(scalefactor); for(size_t i = 0; i < _bins.size(); ++i) _bins[i].scaleW(scalefactor); } + //@} + /// @name Operators //@{ @@ -258,16 +302,18 @@ return _binHashSparse == other._binHashSparse; } + /// Check if the binning of two of the Axis is different bool operator != (const Axis1D& other) const { return ! operator == (other); } + /// Add two Axis together Axis1D<BIN1D,DBN>& operator += (const Axis1D<BIN1D,DBN>& toAdd) { - if(*this != toAdd) throw LogicError("YODA::Histo1D: Cannot add axes with differen binnings."); - - for(size_t i = 0; i < _bins.size(); ++i) { + if (*this != toAdd) throw LogicError("YODA::Histo1D: Cannot add axes with differen binnings."); + + for (size_t i = 0; i < _bins.size(); ++i) { _bins[i] += toAdd.bins().at(i); } @@ -277,11 +323,12 @@ return *this; } - /// Substract two Axis + + /// Subtract two Axis Axis1D<BIN1D,DBN>& operator -= (const Axis1D<BIN1D,DBN>& toSubtract) { - if(*this != toSubtract) throw LogicError("YODA::Histo1D: Cannot add axes with differen binnings."); - - for(size_t i = 0; i < _bins.size(); ++i) { + if (*this != toSubtract) throw LogicError("YODA::Histo1D: Cannot add axes with differen binnings."); + + for (size_t i = 0; i < _bins.size(); ++i) { _bins[i] -= toSubtract.bins().at(i); } @@ -290,55 +337,60 @@ _overflow += toSubtract._overflow; return *this; } - + //@} + private: + /// Makes an axis out of a vector of edges void _mkAxis(const std::vector<double>& binedges) { - for(size_t i = 0; i < binedges.size() - 1; ++i) { - /// If the bin that is to be added doesn't conflict with an - /// existing axis - if(!_findCuts(binedges[i], binedges[i+1])) { - /// Add this bin to the bin vector + for (size_t i = 0; i < binedges.size() - 1; ++i) { + // If the bin that is to be added doesn't conflict with an + // existing axis + if (!_findCuts(binedges[i], binedges[i+1])) { + // Add this bin to the bin vector _bins.push_back(BIN1D(binedges[i], binedges[i+1])); } } - - /// Regenerate the Bin hash + + // Regenerate the Bin hash _mkBinHash(); } + /// Makes an axis from a vector of bins void _mkAxis(const Bins& bins) { _bins = bins; _mkBinHash(); } - + + /// Regenerate the bin hash void _mkBinHash() { - /// Remove the old cache, if it exists + // Remove the old cache, if it exists _binHashSparse.clear(); - /// For each of the bins, add the low and high edge to the cache - for(size_t i = 0; i < _bins.size(); ++i) { + // For each of the bins, add the low and high edge to the cache + for (size_t i = 0; i < _bins.size(); ++i) { _binHashSparse.push_back(std::make_pair(_bins[i].xMin(), i)); _binHashSparse.push_back(std::make_pair(_bins[i].xMax(), i)); } - /// Sort the cache to be nondecreasing to enable - /// fast searches + // Sort the cache to be nondecreasing to enable + // fast searches std::sort(_binHashSparse.begin(), _binHashSparse.end()); } - /// Binary search algorithm. For an in-depth explanation - /// please see the documentation of an analogous one in Axis2D + + // Binary search algorithm. For an in-depth explanation + // please see the documentation of an analogous one in Axis2D size_t _binaryS(double value, size_t lower, size_t higher) const { - if(lower == higher) return lower; + if (lower == higher) return lower; size_t where = (higher+lower)/2; - if(value >= _binHashSparse[where].first) { - if(where == _binHashSparse.size() - 1) return where; + if (value >= _binHashSparse[where].first) { + if (where == _binHashSparse.size() - 1) return where; if(value <= _binHashSparse[where+1].first) return where; return _binaryS(value, where, higher); } @@ -347,29 +399,33 @@ return _binaryS(value, lower, where); } - /// Check if a hypothetical bin to be added, starting at from and + + /// Check if a hypothetical bin to be added, starting at from and /// ending at to will partially overlap any of the existing bins. bool _findCuts(const double& from, const double& to) const { size_t index1 = _binaryS(from, 0, _binHashSparse.size()); size_t index2 = _binaryS(to, 0, _binHashSparse.size()); return !(index1 == index2); } - + + /// Check if there are any gaps in the Axis bool _isGapless(size_t from, size_t to) { size_t start = _binaryS(_bins[from].xMin(), 0, _binHashSparse.size()); size_t end = _binaryS(_bins[to].xMin(), 0, _binHashSparse.size()); - for(size_t i = start; i < end; i++) { + for(size_t i = start; i < end; i++) { if(!fuzzyEquals(_binHashSparse[i].first, _binHashSparse[i+1].first) && _binHashSparse[i].second != _binHashSparse[i+1].second) return false; } - return true; } + + private: + /// @name Data structures //@{ - + /// Bins vector Bins _bins; @@ -377,31 +433,34 @@ DBN _dbn; /// Under- and over- flows - DBN _underflow; DBN _overflow; + DBN _underflow; + DBN _overflow; /// Cached bin edges std::vector<Edge> _binHashSparse; //@} + }; - /// Add the statistics on two axis. + + /// Add the statistics on two axes. template <typename BIN1D, typename DBN> - Axis1D<BIN1D,DBN> operator + (const Axis1D<BIN1D,DBN>& first, const Axis1D<BIN1D,DBN>& second) { + inline Axis1D<BIN1D,DBN> operator + (const Axis1D<BIN1D,DBN>& first, const Axis1D<BIN1D,DBN>& second) { Axis1D<BIN1D,DBN> tmp = first; tmp += second; return tmp; } - + /// Subtract the statistics on two axis. template <typename BIN1D, typename DBN> - Axis1D<BIN1D,DBN> operator - (const Axis1D<BIN1D,DBN>& first, const Axis1D<BIN1D,DBN>& second) { + inline Axis1D<BIN1D,DBN> operator - (const Axis1D<BIN1D,DBN>& first, const Axis1D<BIN1D,DBN>& second) { Axis1D<BIN1D,DBN> tmp = first; tmp -= second; return tmp; } + } #endif - Modified: trunk/include/YODA/Bin1D.h ============================================================================== --- trunk/include/YODA/Bin1D.h Mon Aug 22 18:01:34 2011 (r328) +++ trunk/include/YODA/Bin1D.h Tue Aug 23 01:54:26 2011 (r329) @@ -233,11 +233,23 @@ //@} - protected: - /// @name Named operators //@{ + /// Merge two adjacent bins + Bin1D<DBN>& merge(const Bin1D<DBN>& b) { + if (fuzzyEquals(_edges.second, b._edges.first)) { + _edges.second = b._edges.second; + } else if (fuzzyEquals(_edges.second, b._edges.first)) { + _edges.first = b._edges.first; + } else { + throw LogicError("Attempted to merge two non-adjacent bins"); + } + _dbn += b._dbn; + return *this; + } + + /// Add two bins (internal, explicitly named version) /// /// This operator is defined for adding two bins with equivalent binning. @@ -250,6 +262,7 @@ return *this; } + /// Subtract one bin from another (internal, explicitly named version) /// /// This operator is defined for subtracting two bins with equivalent binning. Modified: trunk/include/YODA/Histo1D.h ============================================================================== --- trunk/include/YODA/Histo1D.h Mon Aug 22 18:01:34 2011 (r328) +++ trunk/include/YODA/Histo1D.h Tue Aug 23 01:54:26 2011 (r329) @@ -132,7 +132,6 @@ } /// Access the bin vector - /// @todo Actually, it's a Histo std::vector<YODA::HistoBin1D>& bins() { return _axis.bins(); } @@ -182,7 +181,7 @@ return _axis.overflow(); } - /// Add a new bin specifying its lower and upper bound + /// Add a new bin specifying its lower and upper bound void addBin(double from, double to) { _axis.addBin(from, to); } @@ -191,7 +190,7 @@ void addBin(std::vector<double> edges) { _axis.addBin(edges); } - + /// Add new bins specifying a beginning and end of each of them void addBin(std::vector<std::pair<double,double> > edges) { _axis.addBin(edges); Modified: trunk/src/Histo1D.cc ============================================================================== --- trunk/src/Histo1D.cc Mon Aug 22 18:01:34 2011 (r328) +++ trunk/src/Histo1D.cc Tue Aug 23 01:54:26 2011 (r329) @@ -144,10 +144,11 @@ ey = y * sqrt( sqr(b1.heightErr()/b1.height()) + sqr(b2.heightErr()/b2.height()) ); break; case BINOMIAL: - /// @todo Check that this is correct + /// @todo Check that this is correct -- isn't it a problem that this varies if the same scale + /// factor is applied to the weights on bins 1 and 2? /// @todo I think this is only valid if the fills of b1 are a subset of the fills of b2. Throw an /// error if Neff(b1) > Neff(b2) - ey = std::sqrt( b1.effNumEntries() * (1- b1.effNumEntries()/b2.effNumEntries()) )/b2.effNumEntries(); + ey = sqrt( b1.effNumEntries() * (1 - b1.effNumEntries()/b2.effNumEntries()) ) / b2.effNumEntries(); break; }
More information about the yoda-svn mailing list |