|
[yoda-svn] r537 - in trunk: . include/YODA pyext pyext/yoda pyext/yoda/includeblackhole at projects.hepforge.org blackhole at projects.hepforge.orgThu Nov 15 16:38:59 GMT 2012
Author: davemallows Date: Thu Nov 15 16:38:59 2012 New Revision: 537 Log: Substantial changes to Axis2D. Axis2D now uses BinSearchers as Axis1D does. Modified: trunk/ChangeLog trunk/include/YODA/Axis2D.h trunk/include/YODA/Bin2D.h trunk/include/YODA/Histo2D.h trunk/include/YODA/Profile2D.h trunk/pyext/setup.py.in trunk/pyext/yoda/axis.pyx trunk/pyext/yoda/core.pyx trunk/pyext/yoda/declarations.pxd trunk/pyext/yoda/include/Axis2D_BIN2D_DBN.pyx trunk/pyext/yoda/include/Bin2D_DBN.pyx trunk/pyext/yoda/include/Histo1D.pyx trunk/pyext/yoda/include/Histo2D.pyx trunk/pyext/yoda/include/Profile2D.pyx Modified: trunk/ChangeLog ============================================================================== --- trunk/ChangeLog Thu Nov 15 16:19:42 2012 (r536) +++ trunk/ChangeLog Thu Nov 15 16:38:59 2012 (r537) @@ -1,3 +1,8 @@ +2012-11-15 Dave Mallows <dave.mallows at gmail.com> + + * Commited numerous changes to Axis2D. Axis2D now uses BinSearcher as + with Axis1D. + 2012-11-15 Andy Buckley <andy.buckley at cern.ch> * Improving division and efficiency treatments, and allowing Modified: trunk/include/YODA/Axis2D.h ============================================================================== --- trunk/include/YODA/Axis2D.h Thu Nov 15 16:19:42 2012 (r536) +++ trunk/include/YODA/Axis2D.h Thu Nov 15 16:38:59 2012 (r537) @@ -1,516 +1,451 @@ -// -*- C++ -*- -// -// This file is part of YODA -- Yet more Objects for Data Analysis -// Copyright (C) 2008-2012 The YODA collaboration (see AUTHORS for details) -// #ifndef YODA_Axis2D_h #define YODA_Axis2D_h #include "YODA/AnalysisObject.h" #include "YODA/Exceptions.h" #include "YODA/Bin.h" -#include "YODA/HistoBin2D.h" -#include "YODA/Utils/cachedvector.h" #include "YODA/Utils/MathUtils.h" -#include "YODA/Dbn2D.h" - -#include <algorithm> +#include "YODA/Utils/BinSearcher.h" #include <limits> -#include <vector> -namespace YODA { +#include <string> +namespace YODA { - /// @brief 1D bin container + /// @brief 2D bin container /// /// This class handles most of the low-level operations on an axis of bins - /// arranged in a 2D grid (including gaps). + /// arranged in a 2D line (including gaps). template <typename BIN2D, typename DBN> class Axis2D { public: - /// @name Typedefs + /// Typedefs //@{ /// Bin type typedef BIN2D Bin; - /// The internal bin container type. Not used for searching. + /// A vector containing 2D bins. Not used for searching. typedef typename std::vector<Bin> Bins; - /// The internal container type for outflow distributions. - typedef std::map<size_t, DBN> Outflows; + // Distinguishing between cuts and pairs is useful + typedef std::vector<double> EdgeCuts; + typedef std::pair<double, double> EdgePair1D; + typedef std::pair<EdgePair1D, EdgePair1D> EdgePair2D; + typedef std::vector<EdgePair2D> EdgePair2Ds; - /// @brief Type used to implement a search table of low bin edges (in 2D) mapped to bin indices. - /// An index of -1 indicates a gap interval, without a corresponding bin. - typedef std::map<double, long int> SubBinHash; - typedef std::map<double, SubBinHash> BinHash; + // Ordered in some arbitrary way: see outflow(int, int) + typedef std::vector<DBN> Outflows; //@} - - /// @name Constructors: + /// @name Constructors //@{ - /// Empty constructor + // Empty constructor Axis2D() - : _isPerfectGrid(true), _locked(false) + : _locked(false) { reset(); } - - /// A constructor with specified x and y axis bin limits. - Axis2D(const std::vector<double>& xedges, const std::vector<double>& yedges) - : _isPerfectGrid(true), _locked(false) + /// A constructor with specified x and y axis bin cuts. + Axis2D(const EdgeCuts &xedges, + const EdgeCuts &yedges) + : _locked(false) { - _addBins(xedges, yedges); + addBins(xedges, yedges); reset(); } - - /// Most standard constructor accepting X/Y ranges and number of bins + /// Constructor accepting X/Y ranges and number of bins /// on each of the axis. Both axes are divided linearly. Axis2D(size_t nbinsX, const std::pair<double,double>& rangeX, size_t nbinsY, const std::pair<double,double>& rangeY) - : _isPerfectGrid(true), _locked(false) + : _locked(false) { - _addBins(linspace(nbinsX, rangeX.first, rangeX.second), - linspace(nbinsY, rangeY.first, rangeY.second)); + addBins(linspace(nbinsX, rangeX.first, rangeX.second), + linspace(nbinsY, rangeY.first, rangeY.second)); reset(); } - /// Constructor accepting a list of bins Axis2D(const Bins& bins) { - _addBins(bins); + addBins(bins); reset(); } - /// State-setting constructor for persistency Axis2D(const Bins& bins, const DBN& totalDbn, const Outflows& outflows) - : _bins(bins), _dbn(totalDbn), _outflows(outflows), - _isPerfectGrid(true), _locked(false) + : _dbn(totalDbn), _outflows(outflows), + _locked(false) // Does this make sense? { if (_outflows.size() != 8) { throw Exception("Axis2D outflow containers must have exactly 8 elements"); } - _updateAxis(); + addBins(bins); } - //@} + void reset() { + _dbn.reset(); + _outflows.assign(8, DBN()); + } + + /// Get the number of bins. + size_t numBins() const { + return _bins.size(); + } + + /// Get the number of bins on the x-axis. This is only sensible for + /// perfectly regular gridded bins. For irregular binnings, this is + /// the number of cuts that were necessary to grid the data. + size_t numBinsX() const { + return _nx; + } + /// Get the number of bins on the y-axis. This is only sensible for + /// perfectly regular gridded bins. For irregular binnings, this is + /// the number of cuts that were necessary to grid the data. + size_t numBinsY() const { + return _ny; + } - /// @name Bin insertion operators + //@} + // + /// @name Statistics accessor functions //@{ - /// @brief Bin addition operator - /// - /// This operator is supplied with the extremal coordinates of just a new - /// bin, which is then constructed and added as usual. - // void addBin(double lowX, double lowY, double highX, double highY) { - void addBin(double, double, double, double) { - /// @todo TODO + /// Get the outflow by x-index and y-index -- e.g. (+1, -1) is outside by + /// being greater than the greatest x-edge and less than the lowest y-edge. + DBN& outflow(int ix, int iy) { + // Lookup table for mapping. This is necessary as there is no + // numerical way to skip the eighth item. This also allows for + // arbitrary orderings. + const static unsigned char outflowMapping[9] = {0, 1, 2, 3, 8, 4, 5, 6, 7}; + ++ix; + ++iy; + if (ix > 2 || iy > 2) + throw UserError( + "Outflow index out of range: valid indices are -1, 0, 1"); + // Find the real index + size_t realindex = outflowMapping[3 * ix + iy]; + // Check we're not using the invalid index + if (realindex == 8) { + throw UserError("(0, 0) is not a valid outflow index"); + } + return _outflows[realindex]; + } - /// @todo Check for overlaps and whether axis is locked + /// Get the outflow by x-index and y-index -- e.g. (+1, -1) is outside by + /// being greater than the greatest x-edge and less than the lowest y-edge. + /// (const version) + const DBN& outflow(int ix, int iy) const { + const static unsigned char outflowMapping[9] = {0, 1, 2, 3, 8, 4, 5, 6, 7}; + ++ix; + ++iy; + if (ix > 2 || iy > 2) + throw UserError( + "Outflow index out of range: valid indices are -1, 0, 1"); + // Find the real index + size_t realindex = outflowMapping[3 * ix + iy]; + // Check we're not using the invalid index + if (realindex == 8) { + throw UserError("(0, 0) is not a valid outflow index"); + } + return _outflows[realindex]; } - //@} + /// Scale each bin as if the entire x-axis had been scaled by this + /// factor. + void scaleX(double xscale) { + scaleXY(xscale, 1.0); + } + + /// Scale each bin as if the entire y-axis had been scaled by this + /// factor. + void scaleY(double yscale) { + scaleXY(1.0, yscale); + } + + /// Scale each bin as if the entire x and y-axes had been scaled by + /// their respective factors. + void scaleXY(double sx, double sy) { + _dbn.scaleXY(sx, sy); + foreach (DBN &dbn, _outflows) + dbn.scaleXY(sx, sy); + foreach (Bin &bin, _bins) + bin.scaleXY(sx, sy); + _updateAxis(_bins); + } - /// @name Modifiers - //@{ + /// Rescale as if all fill weights had been different by factor @a + /// scalefactor. + void scaleW(double scaleFactor) { + _dbn.scaleW(scaleFactor); + foreach (DBN &dbn, _outflows) + dbn.scaleW(scaleFactor); + foreach (Bin &bin, _bins) + bin.scaleW(scaleFactor); + _updateAxis(_bins); + } - // /// Merge a range of bins, given the bin IDs of points at the lower-left and upper-right. - // /// @todo TODO - // void mergeBins(size_t from, size_t to) { - // /// Acquire the starting/ending bins - // BIN2D& start = bin(from); - // BIN2D& end = bin(to); - - // /// Set the bin to be added as a starting bin - // /// and then remove the unneeded starting bin from - // /// the list of bins. - // BIN2D temp = start; - // eraseBin(from); - - // // Sanity-check of input indices - // if (start.midpoint().first > end.midpoint().first) { - // throw RangeError("The start bin has a greater x value than the end bin."); - // } - // if (start.midpoint().second > end.midpoint().second) { - // std::cout << "Start: " << start.midpoint().second; - // std::cout << " end: " << end.midpoint().second << std::endl; - // throw RangeError("The start bin has a greater y value than the end bin."); - // } - - // /// Create a vector that will contain indexes of bins that - // /// will be removed after merging them with our 'main' bin. - // std::vector<size_t> toRemove; - - // /// Look for lower/upper limit of the merge function operation. - // /// i.e.: search for index of lowEdgeY of starting bin in _binHashSparse - // /// and lowEdgeY of ending bin. This way we don't have to scroll through all - // /// the bins to check which ones we have to add. - // for (size_t y = (*_binHashSparse.first._cache.lower_bound(start.yMin())).second; - // y <= (*_binHashSparse.first._cache.lower_bound(end.yMin())).second; ++y) { - // /// For all the bins that are in the bounds specified in previous for - // for (size_t x = 0; x < _binHashSparse.first[y].second.size(); ++x) { - // /// If the bin lies in a rectangle produced by starting and ending bins - // /// (ie. the one spanned by lower-left point of starting bin and top-right - // /// point of ending bin) and was not merged already: - // if (((_binHashSparse.first[y].second[x].second.first > start.xMin() || - // fuzzyEquals(_binHashSparse.first[y].second[x].second.first, start.xMin())) && - // (_binHashSparse.first[y].second[x].second.second < end.xMax() || - // fuzzyEquals(_binHashSparse.first[y].second[x].second.second, end.xMax())))&& - // !(std::find(toRemove.begin(), toRemove.end(), _binHashSparse.first[y].second[x].first) != toRemove.end())) - // { - // /// Merge it with the temp bin and mark it as ready for removal - // temp += bin(_binHashSparse.first[y].second[x].first); - // toRemove.push_back(_binHashSparse.first[y].second[x].first); - // } - // } - // } - - // /// Now, drop the bins to be dropped - // /// Keeping in mind that the bins must be removed from the highest index - // /// down, otherwise we will end up removing other bins that we intend to - // std::sort(toRemove.begin(), toRemove.end()); - // std::reverse(toRemove.begin(), toRemove.end()); - // foreach(size_t remove, toRemove) eraseBin(remove); - - // /// Add edges of our merged bin to _binHashSparse and don't create a default - // /// empty bin. - // _addEdge(temp.edges(), _binHashSparse, false); - - // /// Add the actual merged bin to the Axis. - // _bins.push_back(temp); - - - - // /// And regenerate the caches on _binHashSparse - // _binHashSparse.first.regenCache(); - // _binHashSparse.second.regenCache(); - // } - - // /// Merge a range of bins giving start and end coordinates - // void mergeBins(double startX, double startY, double endX, double endY) { - // mergeBins(binByCoord(startX, startY), binByCoord(endX, endY)); - // } - - // /// Rebin by an interger factor - // /// NOT YET WORKING!!11!!1111 - // void rebin(size_t factorX, size_t factorY) { - // if (!isGrid()) throw GridError("Rebinning by a factor can only act on full grids!"); - // if(factorX < 1 || factorY < 1) throw RangeError("Factors cannot be smaller than unity!"); - - // size_t binsInColumn = _binHashSparse.first.size() - 1; - - // std::cout << std::endl << "Bins in column: " << binsInColumn << std::endl; - // std::cout << "Number of bins: " << _bins.size() << std::endl; - // size_t startIndex = 0; - // while(true) { - // size_t endIndex = startIndex; - // for(size_t i = 1; i < factorY; ++i){ - // if(_hasAbove(endIndex) == 1) endIndex++; - // else break; - // } - // binsInColumn -= endIndex - startIndex; - // for(size_t i = 1; i < factorX; ++i){ - // if(endIndex + binsInColumn < _bins.size()) endIndex += binsInColumn; - // else break; - // } - // if(endIndex + 1 >= _bins.size()) break; - // mergeBins(startIndex, endIndex); - // if(startIndex + 1 < _bins.size()) startIndex++; - // else break; - // if(_hasAbove(startIndex-1) == 0) binsInColumn = _binHashSparse.first.size() -1; - // } + /// Remove the bin at the given index. If many bins need to be + /// removed, prefer eraseBins(vector[size_t] &) over many calls to this, + /// as recreating the binhash is comparatively expensive. + void eraseBin(size_t i) { + if (i >= numBins()) + throw RangeError("Bin index is out of range"); + + // Temporarily unlock the axis during the update + _bins.erase(_bins.begin() + i); + _updateAxis(_bins); + } - // } + /// Erase a rectangle of bins. + void eraseBins(const size_t from, const size_t to) + { + if (from >= numBins()) + throw RangeError("Initial bin index is out of range"); + if (from >= numBins()) + throw RangeError("Final bin index is out of range"); + + Bin &fromBin = bin(from); + Bin &toBin = bin(to); + + eraseBins( + std::make_pair(fromBin.xMin(), toBin.xMax()), + std::make_pair(fromBin.yMin(), toBin.yMax())); + } + + /// Erase bins in an x- and y-range. Any bins which lie entirely within the + /// region are deleted. If any part of the bin lies outside this + /// range, the bin remains, so this has similar behaviour to select + /// tools in vector graphics gui packages. + + // todo: any ideas on how to test this? + void eraseBins(const std::pair<double, double> xrange, + const std::pair<double, double> yrange) + { + size_t xiLow = _binSearcherX.index(xrange.first) - 1; + size_t xiHigh = _binSearcherX.index(xrange.second) - 1; + size_t yiLow = _binSearcherY.index(yrange.first) - 1; + size_t yiHigh = _binSearcherY.index(yrange.second) - 1; - /// Reset the axis statistics - void reset() { - _dbn.reset(); - for (int ix = -1; ix <= 1; ++ix) { - for (int iy = -1; iy <= 1; ++iy) { - if (ix == 0 && iy == 0) continue; - outflow(ix, iy).reset(); + std::vector<bool> deleteMask(numBins(), false); + + for (size_t yi = yiLow; yi < yiHigh; yi++) { + for (size_t xi = xiLow; xi < xiHigh; xi++) { + ssize_t i = _indexes[_index(_nx, xi, yi)]; + if (i == -1 || deleteMask[i]) continue; + if (bin(i).fitsInside(xrange, yrange)) + deleteMask[i] = true; } } - foreach(Bin& bin, _bins) bin.reset(); - _locked = false; + + // Now we just update + eraseBins(deleteMask); } + /// Erase using a vector<bool>, where true represents that a bin + /// will be deleted, and false means it will be kept. + void eraseBins(const std::vector<bool> deleteMask) { + Bins newBins; + for (size_t i = 0; i < numBins; i++) + if (!deleteMask[i]) + newBins.push_back(bins(i)); - /// Set the axis lock state - void _setLock(bool locked) { - _locked = locked; + _update(newBins); } - //@} + //@todo + bool _gapInRange(size_t from, size_t to) { + Bin &toBin = bin(to); + Bin &fromBin = bin(from); + return true; + } - /// @name Accessors - //@{ + //@todo + void rebin(size_t n) { + } - /// Get the total number of bins - /// @todo Can't this just be numBins? - size_t numBins() const { - return _bins.size(); + /// Set the axis lock state + void _setLock(bool locked) { + _locked = locked; } - /// Get the number of bins along X axis - size_t numBinsX() const { - return (numBins() > 0) ? _binhash.size() : 0; + /// Return the lowest-valued bin edge along the x-axis + const double lowEdgeX() const { + return _xRange.first; } - /// Get the number of bins along Y axis - size_t numBinsY() const { - return (numBins() > 0) ? _binhash.begin()->second.size() : 0; + /// Alias for lowEdgeX() + const double xMin() const { + return lowEdgeX(); } + /// Return the highest-valued bin edge along the x-axis + const double highEdgeX() const { + return _xRange.second; + } - /// Get the value of the lowest x-edge on the axis - double lowEdgeX() const { - if (numBins() == 0) throw RangeError("This axis contains no bins and so has no defined range"); - return _binhash.begin()->first; + /// Alias for highEdgeX() + const double xMax() const { + return highEdgeX(); } - /// A alias for lowEdgeX() - double xMin() const { return lowEdgeX();} - /// Get the value of the highest x-edge on the axis - double highEdgeX() const { - if (numBins() == 0) throw RangeError("This axis contains no bins and so has no defined range"); - return (--_binhash.end())->first; + /// Return the lowest-valued bin edge along the y-axis + const double lowEdgeY() const { + return _yRange.first; } - /// Alias for highEdgeX() - double xMax() const { return highEdgeX();} - /// Get the value of the lowest y-edge on the axis - double lowEdgeY() const { - if (numBins() == 0) throw RangeError("This axis contains no bins and so has no defined range"); - return _binhash.begin()->second.begin()->first; - } - /// A alias for lowEdgeY() - double yMin() const { return lowEdgeY();} - - /// Get the value of the highest y-edge on the axis - double highEdgeY() const { - if (numBins() == 0) throw RangeError("This axis contains no bins and so has no defined range"); - return (--_binhash.begin()->second.end())->first; + /// Alias for lowEdgeY() + const double yMin() const { + return lowEdgeY(); } - /// Alias for highEdgeY() - double yMax() const { return highEdgeY();} + /// Return the highest-valued bin edge along the y-axis + const double highEdgeY() const { + return _yRange.second; + } - /// Get the bins (non-const version) - Bins& bins() { - return _bins; + /// Alias for highEdgeY() + const double yMax() const { + return highEdgeY(); } - /// Get the bins (const version) - const Bins& bins() const { - return _bins; + + /// Add a bin, providing its x- and y- edge ranges + void addBin(EdgePair1D xrange, EdgePair1D yrange) { + _checkUnlocked(); + Bins newBins(_bins); + newBins.push_back(Bin(xrange, yrange)); + _updateAxis(newBins); } + /// Add a vector of pre-made bins + void addBins(const Bins &bins) { + if (bins.size() == 0) return; - /// @brief Get an outflow (non-const version) - /// - /// Two indices are used, for x and y: -1 = underflow, 0 = in-range, and +1 = overflow. - /// (0,0) is not a valid overflow index pair, since it is in range for both x and y. - DBN& outflow(size_t ix, size_t iy) { - if (ix == 0 && iy == 0) throw UserError("(0,0) is not a valid Axis2D overflow index"); - if (abs(ix) > 1 || abs(iy) > 1) throw UserError("Axis2D overflow indices are -1, 0, 1"); - size_t realindex = 3*(ix+1) + (iy+1); - return _outflows[realindex]; + _checkUnlocked(); + + Bins newBins(_bins); + + foreach(const Bin &b, bins) { + newBins.push_back(b); + } + + _updateAxis(newBins); } - /// @brief Get an outflow (const version) - /// - /// Two indices are used, for x and y: -1 = underflow, 0 = in-range, and +1 = overflow. - /// (0,0) is not a valid overflow index pair, since it is in range for both x and y. - const DBN& outflow(size_t ix, size_t iy) const { - if (ix == 0 && iy == 0) throw UserError("(0,0) is not a valid Axis2D overflow index"); - if (abs(ix) > 1 || abs(iy) > 1) throw UserError("Axis2D overflow indices are -1, 0, 1"); - size_t realindex = 3*(ix+1) + (iy+1); - return _outflows.find(realindex)->second; + + /// Add a contiguous set of bins to an axis, via their list of edges + void addBins(const std::vector<double> &xcuts, + const std::vector<double> &ycuts) { + if (xcuts.size() == 0) return; + if (ycuts.size() == 0) return; + + _checkUnlocked(); + + Bins newBins(_bins); + + for (size_t xi=0; xi < xcuts.size() - 1; xi++) { + for (size_t yi=0; yi < ycuts.size() - 1; yi++) { + newBins.push_back(Bin( + std::make_pair(xcuts[xi], xcuts[xi+1]), + std::make_pair(ycuts[yi], ycuts[yi+1]))); + } + } + + _updateAxis(newBins); } - /// Get the bin with a given index (non-const version) - BIN2D& bin(size_t index) { - if (index >= _bins.size()) throw RangeError("Bin index out of range."); - return _bins[index]; + /// Get the bin index of the bin containing point (x, y). + ssize_t getBinIndex(double x, double y) const { + size_t xi = _binSearcherX.index(x) - 1; + size_t yi = _binSearcherY.index(y) - 1; + if (xi > _nx) return -1; + if (yi > _ny) return -1; + + return _indexes[_index(_nx, xi, yi)]; } - /// Get the bin with a given index (const version) - const BIN2D& bin(size_t index) const { - if (index >= _bins.size()) throw RangeError("Bin index out of range."); - return _bins[index]; + /// Access bin by index + Bin& bin(size_t i) { + return _bins[i]; } + /// Access bin by index (const) + const Bin& bin(size_t i) const { + return _bins[i]; + } - /// Get a bin at given coordinates (non-const version) - BIN2D& binByCoord(double x, double y) { + /// Get the bin containing point (x, y). + Bin& binByCoord(double x, double y) { const int ret = getBinIndex(x, y); if (ret == -1) throw RangeError("No bin found!!"); return bin(ret); } - /// Get a bin at given coordinates (const version) - const BIN2D& binByCoord(double x, double y) const { + /// Get the bin containing point (x, y) (const). + const Bin& binByCoord(double x, double y) const { const int ret = getBinIndex(x, y); if (ret == -1) throw RangeError("No bin found!!"); return bin(ret); } - - /// Get a bin at given coordinates (non-const version) - BIN2D& binByCoord(std::pair<double, double>& coords) { - return binByCoord(coords.first, coords.second); - } - - /// Get a bin at given coordinates (const version) - BIN2D& binByCoord(std::pair<double, double>& coords) const { - return binByCoord(coords.first, coords.second); - } - - - /// Get the total distribution (non-const version) + /// Return the total distribution (non-const) DBN& totalDbn() { return _dbn; } - /// Get the total distribution (const version) + /// Return the total distribution (const) const DBN& totalDbn() const { return _dbn; } - /// @brief Bin index finder - /// - /// Looks through all the bins to see which one contains the point of - /// interest. - long getBinIndex(double coordX, double coordY) const { - // First check that we are within the axis bounds at all - if (coordX < lowEdgeX() || coordX >= highEdgeX()) return -1; - if (coordY < lowEdgeY() || coordY >= highEdgeY()) return -1; - // Then return the lower-edge lookup (in both directions) from the hash map. - // NB. both upper_bound and lower_bound return values *greater* than (or equal) to coord, - // so we have to step back one iteration to get to the lower-or-equal containing bin edge. - BinHash::const_iterator itaboveX = _binhash.upper_bound(coordX); - SubBinHash::const_iterator itaboveY = (--itaboveX)->second.upper_bound(coordY); - int index = (--itaboveY)->second; - return index; - } - - - // /// Fast column number searcher - // size_t getBinColumn(size_t index) const { - // // Check if assumptions are reasonable - // if (!_isGrid) throw GridError("This operation can only be performed when an array is a grid!"); - // if (index >= _bins.size()) throw RangeError("Index is bigger than the size of bins vector!"); - - // // Find the column number - // size_t ret = (*_binHashSparse.first._cache.lower_bound(approx(bin(index).yMin()))).second; - // return ret; - // } - - - // /// Fast row number searcher - // size_t getBinRow(size_t index) const { - // // Check if assumptions are reasonable - // if (!_isGrid) throw GridError("This operation can only be performed when an array is a grid!"); - // if (index >= _bins.size()) throw RangeError("Index is bigger than the size of bins vector!"); - - // // Find the row number - // size_t ret = (*_binHashSparse.second._cache.lower_bound(approx(bin(index).xMin()))).second; - // return ret; - // } - - - /// @brief Bin eraser - /// - /// Removes a bin at a position. - void eraseBin(size_t index) { - // Check the correctness of assumptions - if (index >= numBins()) throw RangeError("Index is bigger than the size of bins vector!"); - - /// @todo Check for overlaps - /// @todo If within the current bounds, create and rehash - /// @todo If outside the bounds, check for grid compatibility, create and rehash - } - //@} - - - /// @name Scaling operations - //@{ - - /// @brief Rescale the axes - /// - /// Scales the axis with a given scale. - /// @todo Add a specific check for a scaling of 1.0, to avoid doing work when no scaling is wanted. - void scaleXY(double scaleX, double scaleY) { - foreach (Bin bin, _bins) { - bin.scaleXY(scaleX, scaleY); - } - _dbn.scaleXY(scaleX, scaleY); - // Outflows - for (int ix = -1; ix <= 1; ++ix) { - for (int iy = -1; iy <= 1; ++iy) { - if (ix == 0 && iy == 0) continue; - outflow(ix, iy).scaleXY(scaleX, scaleY); - } - } - // Rehash - _updateAxis(); + /// Return the bins vector (non-const) + Bins& bins() { + return _bins; } - - /// Scales the bin weights - void scaleW(double scalefactor) { - foreach (Bin bin, _bins) { - bin.scaleW(scalefactor); - } - _dbn.scaleW(scalefactor); - // Outflows - for (int ix = -1; ix <= 1; ++ix) { - for (int iy = -1; iy <= 1; ++iy) { - if (ix == 0 && iy == 0) continue; - outflow(ix, iy).scaleW(scalefactor); - } - } + /// Return the bins vector (const) + const Bins& bins() const { + return _bins; } - //@} - + /// Equality operator (on binning only) - /// @name Operators - //@{ + // (DM: Doesn't this break the semantics of equality? As it's used only + // rarely, isn't there a real case for having a "binningsCompatible" or + // similar method?) + + bool operator == (const Axis2D& other) const { + for(size_t i=0; i < numBins(); i++) + if (fuzzyEquals(bin(i).lowEdgeX(), other.bin(i).lowEdgeX()) && + fuzzyEquals(bin(i).highEdgeX(), other.bin(i).highEdgeX()) && + fuzzyEquals(bin(i).lowEdgeY(), other.bin(i).lowEdgeY()) && + fuzzyEquals(bin(i).highEdgeY(), other.bin(i).highEdgeY())); + return false; - /// Equality operator (on binning only) - bool operator == (const Axis2D& ) const { - /// @todo TODO return true; } - /// Non-equality operator bool operator != (const Axis2D& other) const { return ! operator == (other); } - /// @brief Addition operator - /// At this stage it is only possible to add histograms with the same binnings. - /// @todo Compatible but not equal binning to come soon. Axis2D<BIN2D, DBN>& operator += (const Axis2D<BIN2D, DBN>& toAdd) { if (*this != toAdd) { throw LogicError("YODA::Axis2D: Cannot add axes with different binnings."); @@ -524,9 +459,6 @@ /// Subtraction operator - /// - /// At this stage it is only possible to subtract histograms with the same binnings. - /// @todo Compatible but not equal binning to come soon. Axis2D<BIN2D, DBN>& operator -= (const Axis2D<BIN2D, DBN>& toSubtract) { if (*this != toSubtract) { throw LogicError("YODA::Axis2D: Cannot add axes with different binnings."); @@ -538,160 +470,123 @@ return *this; } - //@} - - private: - /// Add new bins, constructed from two sorted vectors of edges, to the axis - void _addBins(const std::vector<double>& xbinedges, const std::vector<double>& ybinedges) { - if (_locked) { - throw LockError("Attempting to add bins to a locked axis"); - } - - /// @todo Check that vectors are sorted? - - /// @todo Check whether there is overlap with any existing edges in either direction - // if (_edgeInRange(binedges.front(), binedges.back())) { - // throw RangeError("New bin range overlaps with existing bin edges"); - // } - - // Create and add bins - for (size_t i = 0; i < xbinedges.size() - 1; ++i) { - for (size_t j = 0; j < ybinedges.size() - 1; ++j) { - _bins.push_back( BIN2D(std::make_pair(xbinedges[i], xbinedges[i+1]), - std::make_pair(ybinedges[j], ybinedges[j+1])) ); - } - } - _updateAxis(); - } - - - /// Add new bins to the axis - void _addBins(const Bins& bins) { - if (_locked) { - throw LockError("Attempting to add bins to a locked axis"); - } - for (size_t i = 0; i < bins.size(); ++i) { - - /// @todo Check for 2D edge overlaps - - // if (_edgeInRange(bins[i].xMin(), bins[i].xMax())) { - // throw RangeError("New bin range overlaps with existing bin edges"); - // } - _bins.push_back(bins[i]); - } - _updateAxis(); - } - - - /// Sort the bins vector, and regenerate the bin hash - /// - /// The bin hash is purely for searching, and is generated from the bins list only. - void _updateAxis() { - /// First, set up the collection of low edges based on all unique edges - std::vector<double> xedges, yedges; - std::sort(_bins.begin(), _bins.end()); - for (size_t i = 0; i < numBins(); ++i) { - xedges.push_back(bin(i).xMin()); - xedges.push_back(bin(i).xMax()); // only the unique max edges will "survive" - yedges.push_back(bin(i).yMin()); - yedges.push_back(bin(i).xMax()); // only the unique max edges will "survive" - } - std::sort(xedges.begin(), xedges.end()); - std::vector<double>::iterator itx = std::unique(xedges.begin(), xedges.end()); - xedges.resize(itx - xedges.begin()); - std::sort(yedges.begin(), yedges.end()); - std::vector<double>::iterator ity = std::unique(yedges.begin(), yedges.end()); - yedges.resize(ity - yedges.begin()); - - // Guess that it's a perfect grid if (nxedge-1)*(nyedge-1) == nbins - /// @todo This is not an ideal check: it could be a numerical coincidence. Fix! - _isPerfectGrid = ( (xedges.size()-1)*(yedges.size()-1) == numBins()); - - // Create a double-map hash based on the two sets of low edges. Initialize with null bin indices. - _binhash.clear(); - for (size_t ix = 0; ix < xedges.size() - 1; ++ix) { - _binhash[xedges[ix]] = SubBinHash(); - for (size_t iy = 0; iy < yedges.size() - 1; ++iy) { - _binhash[xedges[ix]][yedges[iy]] = -1; - } - } - - // Loop over bins, setting each one's non-null bin index appropriately in the double-hashmap. - for (size_t ib = 0; ib < numBins(); ++ib) { - - // Find the axis low edges contained within this bin - const double xmin(bin(ib).xMin()), ymin(bin(ib).xMin()); - std::vector<double> xlowedges_in_bin, ylowedges_in_bin; - /// @todo Use std::upper/lower_bound? - for (size_t ix = 0; ix < xedges.size() - 1; ++ix) { - if (xedges[ix] >= xmin) xlowedges_in_bin.push_back(xedges[ix]); - } - /// @todo Use std::upper/lower_bound? - for (size_t iy = 0; iy < yedges.size() - 1; ++iy) { - if (yedges[iy] >= ymin) ylowedges_in_bin.push_back(yedges[iy]); - } - - // Set bin indices - for (std::vector<double>::const_iterator x = xlowedges_in_bin.begin(); x != xlowedges_in_bin.end(); ++x) { - for (std::vector<double>::const_iterator y = ylowedges_in_bin.begin(); y != ylowedges_in_bin.end(); ++y) { - _binhash[*x][*y] = ib; + void _checkUnlocked(void) { + // Ensure that axis is not locked + if (_locked) + throw LockError("Attempting to update a locked axis"); + } + + void _updateAxis(Bins &bins) { + + if (bins.size() == 0) { + _binSearcherX = Utils::BinSearcher(); + _binSearcherY = Utils::BinSearcher(); + _nx = 0; + _ny = 0; + _xRange = std::make_pair(0, 0); + _yRange = std::make_pair(0, 0); + } + + // Sort the bins + std::sort(bins.begin(), bins.end()); + + // Create the cuts + std::vector<double> xcuts, ycuts; + foreach (Bin &bin, bins) { + xcuts.push_back(bin.xMin()); + xcuts.push_back(bin.xMax()); + ycuts.push_back(bin.yMin()); + ycuts.push_back(bin.yMax()); + } + + // Sort the cuts + std::sort(xcuts.begin(), xcuts.end()); + std::sort(ycuts.begin(), ycuts.end()); + + // Get unique elements in x- and y-cuts + // @todo -- fuzzy equality + xcuts.resize(std::unique(xcuts.begin(), xcuts.end()) - xcuts.begin()); + ycuts.resize(std::unique(ycuts.begin(), ycuts.end()) - ycuts.begin()); + + size_t nx = xcuts.size(); + size_t ny = ycuts.size(); + size_t N = nx * ny; + + // Create a sea of gaps + std::vector<ssize_t> indexes(N, -1); + + // Create two BinSearchers + Utils::BinSearcher xSearcher(xcuts); + Utils::BinSearcher ySearcher(ycuts); + + // Iterate through bins and find out which + for(size_t i=0; i < bins.size(); i++) { + Bin &bin = bins[i]; + size_t xiMin= xSearcher.index(bin.xMin()) - 1; + size_t xiMax= xSearcher.index(bin.xMax()) - 1; + + size_t yiMin = ySearcher.index(bin.yMin()) - 1; + size_t yiMax = ySearcher.index(bin.yMax()) - 1; + + for (size_t xi = xiMin; xi < xiMax; xi++) { + for (size_t yi = yiMin; yi < yiMax; yi++) { + size_t ix = _index(nx, xi, yi); + if (indexes[ix] != -1) + throw RangeError("Bin edges overlap!"); + else + indexes[ix] = i; } } - } - // // DEBUG - // for (size_t ix = 0; ix < xedges.size() - 1; ++ix) { - // for (size_t iy = 0; iy < yedges.size() - 1; ++iy) { - // std::cout << xedges[ix] << " " << yedges[iy] << " " << _binhash[xedges[ix]][yedges[iy]] << std::endl; - // } - // } + // Job's a good'n - let's change our class. + _nx = nx; + _ny = ny; - } + _xRange = std::make_pair(xcuts.front(), xcuts.back()); + _yRange = std::make_pair(ycuts.front(), ycuts.back()); + _indexes = indexes; + _bins = bins; - /// @todo Convert from 1D to 2D - // /// Check if there are any bin edges between values @a from and @a to. - // bool _edgeInRange(double from, double to) const { - // return (--_binhash.upper_bound(from)) != (--_binhash.lower_bound(to)); - // } - - - /// @todo Convert from 1D to 2D - // /// Check if there are any gaps in the axis' binning between bin indices @a from and @a to, inclusive. - // bool _gapInRange(size_t ifrom, size_t ito) const { - // assert(ifrom < numBins() && ito < numBins() && ifrom <= ito); - // if (ifrom == ito) return false; - // BinHash::const_iterator start = (--_binhash.upper_bound(bin(ifrom).midpoint())); - // BinHash::const_iterator end = _binhash.upper_bound(bin(ifrom).midpoint()); - // for (BinHash::const_iterator it = start; it != end; ++it) { - // if (it->second == -1) return true; - // } - // return false; - // } + _binSearcherX = xSearcher; + _binSearcherY = ySearcher; + } private: + static size_t _index(size_t nx, size_t x, size_t y) { + return y * nx + x; + } + /// @name Data structures //@{ - /// Bins contained in this histogram + /// Bins vector Bins _bins; /// Total distribution DBN _dbn; - /// Outflow distributions, indexed clockwise from top - Outflows _outflows; + // Outflows + std::vector<DBN> _outflows; - /// Cached bin edges for searching - BinHash _binhash; - - /// Whether the binning exactly matches the hash edges - bool _isPerfectGrid; + // Binsearcher, for searching bins + Utils::BinSearcher _binSearcherX; + Utils::BinSearcher _binSearcherY; + + std::pair<double, double> _xRange; + std::pair<double, double> _yRange; + + // Mapping from binsearcher indices to bin indices (allowing gaps) + std::vector<ssize_t> _indexes; + + // Necessary for bounds checking and indexing + size_t _nx; + size_t _ny; /// Whether modifying bin edges is permitted bool _locked; @@ -700,30 +595,6 @@ }; - - /// @name Operators - //@{ - - /// Addition operator - template <typename BIN2D, typename DBN> - inline Axis2D<BIN2D, DBN> operator + (const Axis2D<BIN2D, DBN>& first, const Axis2D<BIN2D, DBN>& second) { - Axis2D<BIN2D, DBN> tmp = first; - tmp += second; - return tmp; - } - - - /// Subtraction operator - template <typename BIN2D, typename DBN> - inline Axis2D<BIN2D, DBN> operator - (const Axis2D<BIN2D, DBN>& first, const Axis2D<BIN2D, DBN>& second) { - Axis2D<BIN2D, DBN> tmp = first; - tmp -= second; - return tmp; - } - - //@} - - } #endif Modified: trunk/include/YODA/Bin2D.h ============================================================================== --- trunk/include/YODA/Bin2D.h Thu Nov 15 16:19:42 2012 (r536) +++ trunk/include/YODA/Bin2D.h Thu Nov 15 16:38:59 2012 (r537) @@ -332,7 +332,22 @@ //@{ /// Merge two adjacent bins - // @TODO: We still need to add a merge method + /* + Bin2D<DBN>& merge(const Bin2D<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"); + } + // std::cout << "a " << _dbn.sumW() << std::endl; + _dbn += b._dbn; + // std::cout << "b " << _dbn.sumW() << std::endl; + return *this; + } + */ /// Add two bins (internal, explicitly named version) @@ -360,6 +375,33 @@ return *this; } + /// Test whether this bin would fit inside the given area. + bool fitsInside(std::pair<double, double> xrange, + std::pair<double, double> yrange) const + { + return (lowEdgeX() >= xrange.first && + lowEdgeY() >= yrange.first && + highEdgeX() < xrange.second && + highEdgeY() < yrange.second); + } + + /// Test whether a point lies within the current bin + bool bounds(double x, double y) const { + return (x >= lowEdgeX() && x < highEdgeX() + && y >= lowEdgeY() && y < highEdgeY()); + } + + + + /// Test whether two bins are adjacent and, if so, return how as an integer. + int adjacentTo(const Bin2D<DBN> &b) const { + for(int i=0; i<4; i++) { + if (_edges_equal(b, i, (i+2) % 4)) + return i; + } + return -1; + } + //@} @@ -372,10 +414,46 @@ // Distribution of weighted x (and perhaps y) values DBN _dbn; + std::pair<double, double> _edge_par(int i) const { + if (i % 2) + return xedges(); + else + return yedges(); + } + + double _edge_perp(int i) const { + double output; + + switch (i) { + case 0: output = xMax(); break; + case 1: output = yMax(); break; + case 2: output = xMin(); break; + case 3: output = yMin(); break; + } + + return output; + } + + // Check if common edge. + bool _edges_equal(const Bin2D<DBN> &other, const int i, const int j) const { + return other._edges_equal(_edge_perp(i), _edge_par(i), j); + } + + bool _edges_equal(const double perp, + const std::pair<double, double> par, int j) const { + return ( + fuzzyEquals(perp, _edge_perp(j)) && + fuzzyEquals(par.first, _edge_par(j).first) && + fuzzyEquals(par.second, _edge_par(j).second) + ); + } + }; + + /// Add two bins /// /// This "add" operator is defined for adding two bins with equivalent binning. Modified: trunk/include/YODA/Histo2D.h ============================================================================== --- trunk/include/YODA/Histo2D.h Thu Nov 15 16:19:42 2012 (r536) +++ trunk/include/YODA/Histo2D.h Thu Nov 15 16:38:59 2012 (r537) @@ -137,6 +137,28 @@ _axis.scaleXY(scaleX, scaleY); } + + // /// @brief Bin addition operator + // /// Add a bin to an axis described by its x and y ranges. + void addBin(Axis::EdgePair1D xrange, Axis::EdgePair1D yrange) { + _axis.addBin(xrange, yrange); + } + + + /// @brief Bins addition operator + /// Add multiple bins from edge cuts without resetting + void addBins(const Axis::EdgeCuts &xcuts, const Axis::EdgeCuts &ycuts) { + _axis.addBins(xcuts, ycuts); + } + + + /// @brief Bins addition operator + /// Add multiple bins without resetting + void addBins(const Bins &bins) { + _axis.addBins(bins); + } + + // /// Adding bins /// @todo TODO // void addBin(const std::vector<std::pair<std::pair<double,double>, std::pair<double,double> > > coords) { @@ -290,6 +312,7 @@ /// Two indices are used, for x and y: -1 = underflow, 0 = in-range, and +1 = overflow. /// (0,0) is not a valid overflow index pair, since it is in range for both x and y. Dbn2D& outflow(int ix, int iy) { + std::cout << "Histo2D::outflow\n"; return _axis.outflow(ix, iy); } Modified: trunk/include/YODA/Profile2D.h ============================================================================== --- trunk/include/YODA/Profile2D.h Thu Nov 15 16:19:42 2012 (r536) +++ trunk/include/YODA/Profile2D.h Thu Nov 15 16:38:59 2012 (r537) @@ -145,12 +145,26 @@ // } - /// @todo TODO // /// @brief Bin addition operator - // /// Add a bin to an axis described by its lower-left and higher-right point - // void addBin(double lowX, double lowY, double highX, double highY) { - // _axis.addBin(lowX, lowY, highX, highY); - // } + // /// Add a bin to an axis described by its x and y ranges. + void addBin(Axis::EdgePair1D xrange, Axis::EdgePair1D yrange) { + _axis.addBin(xrange, yrange); + } + + + /// @brief Bins addition operator + /// Add multiple bins from edge cuts without resetting + void addBins(const Axis::EdgeCuts &xcuts, const Axis::EdgeCuts &ycuts) { + _axis.addBins(xcuts, ycuts); + } + + + /// @brief Bins addition operator + /// Add multiple bins without resetting + void addBins(const Bins &bins) { + _axis.addBins(bins); + } + /// @todo TODO // /// @brief Bin addition operator Modified: trunk/pyext/setup.py.in ============================================================================== --- trunk/pyext/setup.py.in Thu Nov 15 16:19:42 2012 (r536) +++ trunk/pyext/setup.py.in Thu Nov 15 16:38:59 2012 (r537) @@ -40,12 +40,14 @@ make_templates('Bin1D_DBN', DBN=('Dbn1D', 'Dbn2D')) make_templates('Bin2D_DBN', DBN=('Dbn2D', 'Dbn3D')) +header_files = glob('../include/YODA/*.h') + glob('../include/YODA/Utils/*.h') + extns = [ext('util'), ext('core', statics=[statics], - depends=glob('yoda/include/*.pyx') + depends=glob('yoda/include/*.pyx') + header_files ), ext('axis', statics=[statics], - depends=glob('yoda/include/Axis*.pyx') + depends=glob('yoda/include/Axis*.pyx') + header_files ), ] Modified: trunk/pyext/yoda/axis.pyx ============================================================================== --- trunk/pyext/yoda/axis.pyx Thu Nov 15 16:19:42 2012 (r536) +++ trunk/pyext/yoda/axis.pyx Thu Nov 15 16:38:59 2012 (r537) @@ -11,7 +11,11 @@ # Pure python imports from itertools import repeat, imap from operator import attrgetter -from yoda import * + +from yoda.core import (Dbn1D, Dbn2D, Dbn3D, + HistoBin1D, HistoBin2D, + ProfileBin1D, ProfileBin2D) include "include/Errors.pyx" include "include/Axis1D.pxi" +include "include/Axis2D.pxi" Modified: trunk/pyext/yoda/core.pyx ============================================================================== --- trunk/pyext/yoda/core.pyx Thu Nov 15 16:19:42 2012 (r536) +++ trunk/pyext/yoda/core.pyx Thu Nov 15 16:38:59 2012 (r537) @@ -36,10 +36,8 @@ include "include/Scatter2D.pyx" -#include "include/Bin2D.pxi" -#include "include/ProfileBin2D.pyx" -#include "include/HistoBin2D.pyx" -#include "include/Histo2D.pyx" -#include "include/Profile2D.pyx" - -#include "include/Axis2D.pxi" +include "include/Bin2D.pxi" +include "include/ProfileBin2D.pyx" +include "include/HistoBin2D.pyx" +include "include/Histo2D.pyx" +include "include/Profile2D.pyx" Modified: trunk/pyext/yoda/declarations.pxd ============================================================================== --- trunk/pyext/yoda/declarations.pxd Thu Nov 15 16:19:42 2012 (r536) +++ trunk/pyext/yoda/declarations.pxd Thu Nov 15 16:38:59 2012 (r537) @@ -348,6 +348,8 @@ Bin2D operator + (Bin2D) except+ err Bin2D operator - (Bin2D) except+ err + int adjacentTo(Bin2D) except+ err + ctypedef Bin2D[Dbn2D] Bin2D_Dbn2D ctypedef Bin2D[Dbn3D] Bin2D_Dbn3D # }}} Bin2D @@ -621,8 +623,8 @@ #void mergeBins(size_t, size_t) except+ err #void rebin(int n) except+ err - # void addBin(double, double) except+ err - # void addBins(vector[double] edges) except+ err + void addBin(pair[double, double], pair[double, double]) except+ err + void addBins(vector[double], vector[double]) except+ err # void eraseBin(size_t index) except+ err size_t numBins() except+ err @@ -725,6 +727,10 @@ HistoBin2D bin(size_t ix) except+ err HistoBin2D binByCoord(double x) except+ err + void addBins(vector[HistoBin2D]&) + + void addBin(pair[double, double], pair[double, double]) + # These must be treated as references or the semantics is wrong. # However, these can also throw exceptions! But the two cannot mix, or # Cython puts out rubbish C++. Since this *is* a reported 'major' bug, @@ -743,6 +749,8 @@ double highEdgeX() except+ err double highEdgeY() except+ err + + int findBinIndex(double, double) size_t numBinsX() size_t numBinsY() @@ -820,3 +828,26 @@ void eraseBin(size_t index) except+ err void mergeBins(size_t, size_t) except+ err # Axis1D }}} + +# Axis2D {{{ +cdef extern from "YODA/Axis2D.h" namespace "YODA": + cdef cppclass Axis2D[BIN2D, DBN]: + Axis2D() except+ err + Axis2D(vector[double], vector[double]) except+ err + Axis2D(size_t, pair[double, double], size_t, pair[double, double]) except+ err + Axis2D(vector[BIN2D] bins) except+ err + void addBin(pair[double, double], pair[double, double]) except+ err + size_t numBins() except+ err + vector[BIN2D] &bins() + double lowEdgeX() except+ err + double highEdgeX() except+ err + double lowEdgeY() except+ err + double highEdgeY() except+ err + long getBinIndex(double, double) + void reset() + DBN &totalDbn() + DBN &outflow(int, int) + void eraseBin(size_t index) except+ err + void mergeBins(size_t, size_t) except+ err +# Axis2D }}} + Modified: trunk/pyext/yoda/include/Axis2D_BIN2D_DBN.pyx ============================================================================== --- trunk/pyext/yoda/include/Axis2D_BIN2D_DBN.pyx Thu Nov 15 16:19:42 2012 (r536) +++ trunk/pyext/yoda/include/Axis2D_BIN2D_DBN.pyx Thu Nov 15 16:38:59 2012 (r537) @@ -24,6 +24,9 @@ return util.new_owned_cls( ${DBN}, new c.${DBN}(self._Axis2D().totalDbn())) + def addBin(self, a, b, c, d): + self._Axis1D().addBin(a, b, c, d) + @property def outflow(self, ix, iy): return util.new_owned_cls( @@ -31,9 +34,9 @@ @property def edges(self): - return XY( - Edge(self._Axis2D().lowEdgeX(), self._Axis2D().highEdgeX()), - Edge(self._Axis2D().lowEdgeY(), self._Axis2D().highEdgeY()) + return util.XY( + util.Edges(self._Axis2D().lowEdgeX(), self._Axis2D().highEdgeX()), + util.Edges(self._Axis2D().lowEdgeY(), self._Axis2D().highEdgeY()) ) def reset(self): Modified: trunk/pyext/yoda/include/Bin2D_DBN.pyx ============================================================================== --- trunk/pyext/yoda/include/Bin2D_DBN.pyx Thu Nov 15 16:19:42 2012 (r536) +++ trunk/pyext/yoda/include/Bin2D_DBN.pyx Thu Nov 15 16:38:59 2012 (r537) @@ -22,40 +22,40 @@ """The lower and upper edges.""" cdef pair[double, double] x = self._Bin2D().xedges() cdef pair[double, double] y = self._Bin2D().yedges() - return XY(util.Edges(x.first, x.second), + return util.XY(util.Edges(x.first, x.second), util.Edges(y.first, y.second)) @property def widths(self): """The widths of this bin in the x- and y-dimensions.""" - return XY(self._Bin2D().widthX(), self._Bin2D().widthY()) + return util.XY(self._Bin2D().widthX(), self._Bin2D().widthY()) @property def focus(self): """The focus of the bin in the x- and y-dimensions""" cdef pair[double, double] f = self._Bin2D().focus() - return XY(f.first, f.second) + return util.XY(f.first, f.second) @property def midpoint(self): cdef pair[double, double] f = self._Bin2D().midpoint() - return XY(f.first, f.second) + return util.XY(f.first, f.second) @property def mean(self): - return XY(self._Bin2D().xMean(), self._Bin2D().yMean()) + return util.XY(self._Bin2D().xMean(), self._Bin2D().yMean()) @property def std_dev(self): - return XY(self._Bin2D().xStdDev(), self._Bin2D().yStdDev()) + return util.XY(self._Bin2D().xStdDev(), self._Bin2D().yStdDev()) @property def std_err(self): - return XY(self._Bin2D().xStdErr(), self._Bin2D().yStdErr()) + return util.XY(self._Bin2D().xStdErr(), self._Bin2D().yStdErr()) @property def rms(self): - return XY(self._Bin2D().xRMS(), self._Bin2D().yRMS()) + return util.XY(self._Bin2D().xRMS(), self._Bin2D().yRMS()) # Raw statistics # ################## @@ -99,6 +99,9 @@ # self._Bin2D().merge(deref(other._Bin2D())) # return self + def adjacentTo(Bin2D_${DBN} self, Bin2D_${DBN} other): + return self._Bin2D().adjacentTo(deref(other._Bin2D())) + def __add__(Bin2D_${DBN} self, Bin2D_${DBN} other): return util.new_owned_cls( Bin2D_${DBN}, Modified: trunk/pyext/yoda/include/Histo1D.pyx ============================================================================== --- trunk/pyext/yoda/include/Histo1D.pyx Thu Nov 15 16:19:42 2012 (r536) +++ trunk/pyext/yoda/include/Histo1D.pyx Thu Nov 15 16:38:59 2012 (r537) @@ -82,6 +82,15 @@ """ self._Histo1D().fill(x, weight) + ## HACK hACK HACK HACK HACK + def fill_many(self, xs, weight=1.0): + """ + (x, weight=1.0) -> None. Fill with given x and optional weight. + + """ + for x in xs: + self._Histo1D().fill(x, weight) + def copy(self, char *path=""): """(path="") -> Histo1D. Clone this Histo1D with optional new path.""" return util.new_owned_cls(Histo1D, Modified: trunk/pyext/yoda/include/Histo2D.pyx ============================================================================== --- trunk/pyext/yoda/include/Histo2D.pyx Thu Nov 15 16:19:42 2012 (r536) +++ trunk/pyext/yoda/include/Histo2D.pyx Thu Nov 15 16:38:59 2012 (r537) @@ -35,9 +35,13 @@ def __repr__(self): return "Histo2D" - def fill(self, x, y, weight=1.0): + def fill(self, double x, double y, weight=1.0): self._Histo2D().fill(x, y, weight) + def fill_many(self, xs, ys, weight=1.0): + for x, y in izip(xs, ys): + self._Histo2D().fill(x, y, weight) + def copy(self, char *path=""): return util.new_owned_cls(Histo2D, new c.Histo2D(deref(self._Histo2D()), string(path))) @@ -58,22 +62,22 @@ return self._Histo2D().sumW(overflows) def mean(self, overflows=True): - return XY( + return util.XY( self._Histo2D().xMean(overflows), self._Histo2D().yMean(overflows)) def variance(self, overflows=True): - return XY( + return util.XY( self._Histo2D().xVariance(overflows), self._Histo2D().yVariance(overflows)) def std_dev(self, overflows=True): - return XY( + return util.XY( self._Histo2D().xStdDev(overflows), self._Histo2D().yStdDev(overflows)) def std_err(self, overflows=True): - return XY( + return util.XY( self._Histo2D().xStdErr(overflows), self._Histo2D().yStdErr(overflows)) @@ -99,11 +103,23 @@ #def rebin(self, int n): # self._Histo2D().rebin(n) - #def addBin(self, double low, double high): - # """Add a bin to the Histo2D""" - # self._Histo2D().addBin(low, high) - # return self + def addBin(self, xlow, xhigh, ylow, yhigh): + """Add a bin to the Histo2D""" + self._Histo2D().addBin(pair[double, double](xlow, xhigh), + pair[double, double](ylow, yhigh)) + return self + + def addBins(self, bounds): + v = new vector[HistoBin2D]() + + for xlow, xhigh, ylow, yhigh in bounds: + v.push_back(HistoBin2D(xlow, xhigh, ylow, yhigh)) + + self._Histo2D().addBins(deref(v)) + del v + + # Need to look at all the possible things here... #def addBins(self, edges): # cdef vector[double] cedges # for i in edges: Modified: trunk/pyext/yoda/include/Profile2D.pyx ============================================================================== --- trunk/pyext/yoda/include/Profile2D.pyx Thu Nov 15 16:19:42 2012 (r536) +++ trunk/pyext/yoda/include/Profile2D.pyx Thu Nov 15 16:38:59 2012 (r537) @@ -30,9 +30,6 @@ return util.new_borrowed_cls( ProfileBin2D, & self._Profile2D().bins().at(i), self) - def __repr__(self): - return "<Profile2D at %x>" % id(self) - def fill(self, double x, double y, double z, double weight=1.0): self._Profile2D().fill(x, y, z, weight) @@ -80,7 +77,6 @@ return util.new_borrowed_cls( Dbn3D, new c.Dbn3D(self._Profile2D().totalDbn()), self) - @property def outflow(self, ix, iy): return util.new_borrowed_cls( Dbn3D, new c.Dbn3D(self._Profile2D().outflow(ix, iy)), self) @@ -135,15 +131,20 @@ #def rebin(self, int n): # self._Profile2D().rebin(n) - #def add_bin(self, double low, double high): - # """Add a bin to the Profile2D""" - # self._Profile2D().addBin(low, high) - # return self - - #def add_bins(self, edges): - # cdef vector[double] cedges - # for i in edges: - # cedges.push_back(i) - # self._Profile2D().addBins(cedges) - # return self + def addBin(self, double xlow, double xhigh, double ylow, double yhigh): + """Add a bin to the Profile2D""" + self._Profile2D().addBin(pair[double, double](xlow, xhigh), + pair[double, double](ylow, yhigh)) + return self + + def addBins(self, xcuts, ycuts): + cdef vector[double] _xcuts + cdef vector[double] _ycuts + for x in xcuts: + _xcuts.push_back(x) + for y in ycuts: + _ycuts.push_back(y) + + self._Profile2D().addBins(_xcuts, _ycuts) + return self
More information about the yoda-svn mailing list |