|
[yoda-svn] r523 - in trunk: . include/YODA include/YODA/Utils pyext pyext/yoda pyext/yoda/includeblackhole at projects.hepforge.org blackhole at projects.hepforge.orgTue Aug 7 02:08:16 BST 2012
Author: davemallows Date: Tue Aug 7 02:08:15 2012 New Revision: 523 Log: Converted Axis1D to use Utils::BinSearcher. BinSearcher is a linear lookup for the common regular (lin/log) and near-regular case, which should be the most common cases. BinSearcher falls back to bisection search in the non-regular cases. Added: trunk/include/YODA/Utils/BinSearcher.h trunk/include/YODA/Utils/fastlog.h trunk/pyext/yoda/axis.pyx Modified: trunk/ChangeLog trunk/include/YODA/Axis1D.h trunk/pyext/Makefile.am trunk/pyext/setup.py.in trunk/pyext/yoda/core.pyx trunk/pyext/yoda/declarations.pxd trunk/pyext/yoda/include/Axis1D_BIN1D_DBN.pyx Modified: trunk/ChangeLog ============================================================================== --- trunk/ChangeLog Mon Aug 6 14:49:55 2012 (r522) +++ trunk/ChangeLog Tue Aug 7 02:08:15 2012 (r523) @@ -1,3 +1,6 @@ +2012-08-07 Dave Mallows <dave.mallows at gmail.com> + * Converted Axis1D to use new Utils/BinSearcher. + 2012-08-02 Dave Mallows <dave.mallows at gmail.com> * Heavily refactored Cython bindings Modified: trunk/include/YODA/Axis1D.h ============================================================================== --- trunk/include/YODA/Axis1D.h Mon Aug 6 14:49:55 2012 (r522) +++ trunk/include/YODA/Axis1D.h Tue Aug 7 02:08:15 2012 (r523) @@ -5,12 +5,13 @@ #include "YODA/Exceptions.h" #include "YODA/Bin.h" #include "YODA/Utils/MathUtils.h" +#include "YODA/Utils/BinSearcher.h" +#include <limits> #include <string> namespace YODA { - /// @brief 1D bin container /// /// This class handles most of the low-level operations on an axis of bins @@ -28,13 +29,8 @@ /// A vector containing 1D bins. Not used for searching. typedef typename std::vector<Bin> Bins; - /// @brief Type used to implement a search table of low bin edges mapped to bin indices. - /// An index of -1 indicates a gap interval, without a corresponding bin. - typedef std::map<double, long int> BinHash; - //@} - /// @name Constructors //@{ @@ -48,7 +44,19 @@ Axis1D(const std::vector<double>& binedges) : _locked(false) { - _addBins(binedges); + addBins(binedges); + } + + + /// @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 + /// the statistics + Axis1D(const std::vector<BIN1D>& bins) + : _locked(false) + { + addBins(bins); } @@ -57,7 +65,7 @@ Axis1D(size_t nbins, double lower, double upper) : _locked(false) { - _addBins(linspace(lower, upper, nbins)); + addBins(linspace(lower, upper, nbins)); } @@ -66,23 +74,12 @@ /// Note that not only dimensions of these bins will be set, /// all the contents of the bins will be copied across, including /// the statistics - Axis1D(const std::vector<BIN1D>& bins) - : _locked(false) - { - _addBins(bins); - } - - - /// State-setting constructor for persistency 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), _locked(false) { - _addBins(bins); + addBins(bins); } - //@} - - /// @name Statistics accessor functions //@{ @@ -91,15 +88,15 @@ return bins().size(); } - /// Return a vector of bins (non-const) - Bins& bins() { - return _bins; - } - /// Return a vector of bins (const) const Bins& bins() const { return _bins; } + + /// Return a vector of bins (non-const) + Bins& bins() { + return _bins; + } /// Return the lowest-value bin edge on the axis double lowEdge() const { @@ -143,19 +140,15 @@ return bin(index); } - /// Return the total distribution (non-const) - DBN& totalDbn() { - return _dbn; - } /// Return the total distribution (const) const DBN& totalDbn() const { return _dbn; } - /// Return underflow (non-const) - DBN& underflow() { - return _underflow; + /// Return the total distribution (non-const) + DBN& totalDbn() { + return _dbn; } /// Return underflow (const) @@ -163,9 +156,9 @@ return _underflow; } - /// Return overflow (non-const) - DBN& overflow() { - return _overflow; + /// Return underflow (non-const) + DBN& underflow() { + return _underflow; } /// Return overflow (const) @@ -173,6 +166,12 @@ return _overflow; } + /// Return overflow (non-const) + DBN& overflow() { + return _overflow; + } + + //@} @@ -180,18 +179,11 @@ //@{ /// Returns an index of a bin at a given coord, -1 if no bin matches - long int getBinIndex(double coord) const { - // First check that we are within the axis bounds at all - if (_binhash.empty() || coord < lowEdge() || coord >= highEdge()) return -1; - // Then return the lower-edge lookup 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 itabove = _binhash.upper_bound(coord); - long int index = (--itabove)->second; - return index; + long getBinIndex(double coord) const { + // Yes, this is robust even with an empty axis. + return _indexes[_binsearcher.index(coord)]; } - /// Reset all the bin statistics on the axis void reset() { _dbn.reset(); @@ -212,73 +204,124 @@ /// Merge a series of bins, between the bins identified by indices @a from and @a to void mergeBins(size_t from, size_t to) { // Correctness checking - if (from >= numBins()) throw RangeError("First index is out of range!"); - if ( to >= numBins()) throw RangeError("Second index is out of range!"); - if (_bins[from].xMin() > _bins[to].xMin()) throw RangeError("The starting bin is greater than ending bin!"); - if (_gapInRange(from, to)) throw RangeError("Bin ranges containing binning gaps cannot be merged!"); + if (from >= numBins()) + throw RangeError("Initial merge index is out of range"); + if (to >= numBins()) + throw RangeError("Final merge index is out of range"); + if (from > to) + throw RangeError("Final bin must be greater than initial bin"); + if (_gapInRange(from, to)) + throw RangeError("Bin ranges containing gaps cannot be merged"); + + Bin &b = bin(from); - BIN1D& b = _bins[from]; - for (size_t i = from+1; i <= to; ++i) { + for (size_t i = from + 1; i <= to; i++) b.merge(_bins[i]); - } + eraseBins(from+1, to); } - /// Merge every group of @a n bins, starting from the LHS - void rebin(int n) { - size_t m = 0; - while (m < numBins()) { + void rebin(size_t n) { + for(size_t m=0; m < numBins(); m++) { const size_t end = (m + n - 1 < numBins()) ? m + n -1 : numBins() - 1; if (end > m) mergeBins(m, end); - m += 1; } } - /// Add a bin, providing its low and high edge void addBin(double low, double high) { - std::vector<double> edges; - edges.push_back(low); - edges.push_back(high); - _addBins(edges); + Bins newBins(_bins); + newBins.push_back(Bin(low, high)); + _updateAxis(newBins); } /// Add a contiguous set of bins to an axis, via their list of edges - void addBins(const std::vector<double>& binedges) { - _addBins(binedges); + void addBins(const std::vector<double> &binedges) { + Bins newBins(_bins); + if (binedges.size() == 0) + return; + + double low = binedges.front(); + + for (size_t i=1; i < binedges.size(); i++) { + double high = binedges[i]; + newBins.push_back(Bin(low, high)); + low = high; + } + + _updateAxis(newBins); } + // (The following proposed method would drastically simplify the + // Python interfaces handling of bin adding -- it would also + // simplify some of the hurdles of bins not having default + // constructors) - /// Remove a bin - void eraseBin(size_t index) { - if (index >= _bins.size()) throw RangeError("Index out of range!"); - _bins.erase(_bins.begin() + index); - _updateAxis(); + /// Add a list of bins as pairs of lowEdge, highEdge + void addBins(const std::vector<std::pair<double, double> > &binpairs) { + // Make a copy of the current binning + Bins newBins(_bins); + + // Iterate over given bins + for (size_t i=0; i<binpairs.size(); i++) { + std::pair<double, double> b = binpairs[i]; + newBins.push_back(Bin(b.first, b.second)); + } + _updateAxis(newBins); + } + + + void addBins(const Bins &bins) { + Bins newBins(_bins); + foreach(const Bin &b, bins) { + newBins.push_back(b); + } + _updateAxis(newBins); } + /// Remove a bin + void eraseBin(const size_t i) { + // Might as well erase from the internal bins, as we can guarantee + // consistency. + if (i >= numBins()) + throw RangeError("Bin index is out of range"); + + const bool wasLocked = _locked; + _locked=false; + _bins.erase(_bins.begin() + i); + _updateAxis(_bins); + _locked = wasLocked; + } /// Remove a bin range - void eraseBins(size_t from, size_t to) { - if (from >= numBins()) throw ("First index is out of range!"); - if ( 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!"); + void eraseBins(const size_t from, const size_t to) { + if (from >= numBins()) + throw RangeError("Initial index out of range"); + if (to >= numBins()) + throw RangeError("Final index out of range"); + if (from > to) + throw RangeError("Final index is less than initial index"); + + const bool wasLocked = _locked; + _locked=false; _bins.erase(_bins.begin() + from, _bins.begin() + to + 1); - _updateAxis(); + _updateAxis(_bins); + _locked = wasLocked; } - /// Scale the size of an axis by a factor + // @todo What if somebody passes in negative scalefactor? (test idea) void scaleX(double scalefactor) { _dbn.scaleX(scalefactor); _underflow.scaleX(scalefactor); _overflow.scaleX(scalefactor); - for (size_t i = 0; i < _bins.size(); ++i) _bins[i].scaleX(scalefactor); - _updateAxis(); + for (size_t i = 0; i < _bins.size(); ++i) + _bins[i].scaleX(scalefactor); + _updateAxis(_bins); } - /// Scale the amount of fills by a factor void scaleW(double scalefactor) { _dbn.scaleW(scalefactor); @@ -289,13 +332,17 @@ //@} - /// @name Operators //@{ /// Check if two of the Axis have the same binning bool operator == (const Axis1D& other) const { - return _binhash == other._binhash; + for(size_t i=0; i < numBins(); i++) + if (bin(i).lowEdge() != other.bin(i).lowEdge() || + bin(i).highEdge() != other.bin(i).highEdge()) + return false; + + return true; } @@ -339,75 +386,69 @@ private: - - /// Add new bins, constructed from a sorted vector of edges, to the axis - void _addBins(const std::vector<double>& binedges) { - //std::sort(binedges.begin(), binedges.end()); - if (_locked) { - throw LockError("Attempting to add bins to a locked axis"); - } - if (_edgeInRange(binedges.front(), binedges.back())) { - throw RangeError("New bin range overlaps with existing bin edges"); - } - for (size_t i = 0; i < binedges.size() - 1; ++i) { - _bins.push_back(BIN1D(binedges[i], binedges[i+1])); - } - _updateAxis(); - } - - - /// Add new bins to the axis - void _addBins(const Bins& bins) { + /// Sort the given bins vector, and regenerate the bin searcher + // + /// The bin searcher is purely for searching, and is generated from + /// the bins list only. + void _updateAxis(Bins &bins) { + // Ensure that axis is not locked if (_locked) { - throw LockError("Attempting to add bins to a locked axis"); - } - for (size_t i = 0; i < bins.size(); ++i) { - if (_edgeInRange(bins[i].xMin(), bins[i].xMax())) { - throw RangeError("New bin range overlaps with existing bin edges"); - } - _bins.push_back(bins[i]); + throw LockError("Attempting to update a locked axis"); } - _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() { - std::sort(_bins.begin(), _bins.end()); - _binhash.clear(); - for (size_t i = 0; i < numBins(); ++i) { - // Add low edge hash for each bin - _binhash[bin(i).xMin()] = i; - // If the next bin is not contiguous, add a gap index for the high edge of this bin - if (i+1 < numBins() && !fuzzyEquals(bin(i).xMax(), bin(i+1).xMin())) { - _binhash[bin(i).xMax()] = -1; + // Define the new cuts and indexes + std::vector<double> edgeCuts; + std::vector<long> indexes; + + // Sort the bins + std::sort(bins.begin(), bins.end()); + + // Keep a lag of the last high edge + double last_high = -std::numeric_limits<double>::infinity(); + + for (size_t i=0; i < bins.size(); i++) { + Bin ¤tBin = bins[i]; + + if(currentBin.lowEdge() < last_high) { + throw RangeError("Bin edges overlap"); + } else if (last_high < currentBin.lowEdge()) { + indexes.push_back(-1); + edgeCuts.push_back(currentBin.lowEdge()); } - } - } - - /// Check if there are any bin edges between values @a from and @a to. - bool _edgeInRange(double from, double to) const { - if ( _binhash.empty() ) return false; - return (--_binhash.upper_bound(from)) != (--_binhash.lower_bound(to)); + // Bins check that they are not zero or negative width. It's perfectly + // okay for them to throw an exception here, as we haven't changed + // anything yet. + indexes.push_back(i); + edgeCuts.push_back(currentBin.highEdge()); + + last_high = currentBin.highEdge(); + } + indexes.push_back(-1); + + // Everything was okay, so let's make our changes + //std::cout << "In" << std::endl; + _binsearcher = Utils::BinSearcher(edgeCuts); + //std::cout << "Out" << std::endl; + _indexes = indexes; + _bins = bins; } /// 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; - } + bool _gapInRange(size_t from, size_t to) const { + assert(from < numBins() && to < numBins() && from < to); + if (from == to) return false; + + const size_t from_ix = _binsearcher.index(bin(from).lowEdge()); + const size_t to_ix = _binsearcher.index(bin(to).lowEdge()); + + for (size_t i = from_ix; i <= to_ix; i++) + if (_indexes[i] == -1) + return true; + return false; } - private: /// @name Data structures @@ -423,8 +464,10 @@ DBN _underflow; DBN _overflow; - /// Cached bin edges for searching - BinHash _binhash; + // Binsearcher, for searching bins + Utils::BinSearcher _binsearcher; + // Mapping from binsearcher indices to bin indices (allowing gaps) + std::vector<long> _indexes; /// Whether modifying bin edges is permitted bool _locked; @@ -450,7 +493,6 @@ return tmp; } - } #endif Added: trunk/include/YODA/Utils/BinSearcher.h ============================================================================== --- /dev/null 00:00:00 1970 (empty, because file is newly added) +++ trunk/include/YODA/Utils/BinSearcher.h Tue Aug 7 02:08:15 2012 (r523) @@ -0,0 +1,248 @@ +#ifndef YODA_BINSEARCHER_H +#define YODA_BINSEARCHER_H + +#include <cstdlib> +#include <vector> +#include <limits> +#include <iostream> +#include "YODA/Utils/fastlog.h" +#include <boost/shared_ptr.hpp> +#include <boost/shared_array.hpp> + +namespace YODA { + namespace Utils { + const size_t SEARCH_SIZE = 16; + const size_t BISECT_LINEAR_THRESHOLD = 32; + + + /// @brief Bin estimator + /// + /// This class handles guessing the right bin for a given value. The better + /// the guess, the less time spent looking. + class Estimator { + public: + virtual const double operator() (double x) const { return x; } + }; + + /// @brief Linear bin estimator + /// + /// This class handles guessing a bin on a linear scale. + class LinEstimator : public Estimator { + public: + LinEstimator(double xlow, double xhigh, size_t N) { + _c = xlow; + _max = N + 1; + _m = (double) N / (xhigh - xlow); + } + + const double operator() (double x) const { + return _m * (x - _c); + } + + protected: + double _c, _m; + size_t _max; + }; + + /// @brief Linear bin estimator + /// + /// This class handles guessing a bin on a logarithmic scale. + class LogEstimator : public Estimator { + public: + LogEstimator(double xlow, double xhigh, size_t N) { + _c = log2(xlow); + _max = N + 1; + _m = N / (log2(xhigh) - _c); + } + + const double operator() (double x) const { + return 1 + _m * (fastlog2(x) - _c); + } + + protected: + double _c, _m; + size_t _max, _half; + }; + + + /// @brief Bin searcher + /// + /// Handles low-level bin lookups using a hybrid algorithm that is + /// considerably faster for regular (logarithmic or linear) and near-regular + /// binnings. Comparable performance for irregular binnings. + + // The reason this works is that linear search is faster than bisection + // search up to about 32-64 elements. So we make a guess, and we then do a + // linear search. If that fails, then we bisect on the remainder, + // terminating once bisection search has got the range down to about 32. So + // we actually pay for the fanciness of predicting the bin out of speeding + // up the bisection search by finishing it with a linear search. So in most + // cases, we get constant-time lookups regardless of the space. + + class BinSearcher { + public: + + BinSearcher() { + _est.reset(new LinEstimator(0, 1, 1)); + } + + // Copy constructor + BinSearcher(const BinSearcher& bs) { + _est = bs._est; + _lows = bs._lows; + _max = bs._max; + } + + // Explicit constructor + BinSearcher(std::vector<double> &lows, bool log) + { + _updateLows(lows); + // Use an array for c-style hackery + // (Templates make no sense here) + if (log) { + _est.reset(new LogEstimator(lows.front(), lows.back(), lows.size() - 1)); + } else { + _est.reset(new LinEstimator(lows.front(), lows.back(), lows.size() - 1)); + } + } + + // Fully automatic constructor: give bins and it does the rest! + BinSearcher(std::vector<double> &lows) + { + _updateLows(lows); + + if (!lows.size()) + _est.reset(new LinEstimator(0, 0, 0)); + else { + LinEstimator linEst(lows.front(), lows.back(), lows.size() - 1); + LogEstimator logEst(lows.front(), lows.back(), lows.size() - 1); + + double logsum = 0; + for (size_t i = 0; i < lows.size(); i++) { + logsum += fabs(logEst(lows[i]) - i); + } + + double linsum = 0; + for (size_t i = 0; i < lows.size(); i++) { + linsum += fabs(linEst(lows[i]) - i); + } + + // Use an array for c-style hackery + // (Templates make no sense here) + double log_avg = logsum / lows.size(); + double lin_avg = linsum / lows.size(); + + // This also implicitly works for NaN returned from the log + // There is a subtle bug here if the if statement is the other way + // around, as (nan < linsum) -> false always. + // But (nan > linsum) -> false also. + if (log_avg < lin_avg) { + _est.reset( + new LogEstimator( + lows.front(), lows.back(), lows.size() - 1)); + } else { + _est.reset( + new LinEstimator( + lows.front(), lows.back(), lows.size() - 1)); + } + } + + } + + // Lookup a bin + __attribute__((noinline)) + size_t const index(double x) const { + + size_t index = estimate(x); + + if (x >= _lows[index]) { + size_t di = _linsearch_forward(_lows.get() + index, x, SEARCH_SIZE); + index += di; + if (di == SEARCH_SIZE) { + index = _bisect(_lows.get(), x, index, _max); + } + } else { + size_t di = _linsearch_backward(_lows.get() + index, x, SEARCH_SIZE); + index -= di; + if (di == SEARCH_SIZE) { + index = _bisect(_lows.get(), x, 0, index + 1); + } + } + + return index - 1; + } + + const size_t estimate(const double x) const { + double y = (*_est)(x); + size_t yi = 0 ? (y < 0) : y; + return (yi > _max) ? _max : yi; + } + + const std::pair<double, double> indexRange(size_t ix) const { + return std::make_pair(_lows[ix], _lows[ix+1]); + } + + const bool inRange (std::pair<double, double> range, double x) const { + return (range.first <= x) && (x < range.second); + } + + protected: + void _updateLows(std::vector<double> &lows) { + _lows.reset(new double[lows.size() + 2]); + + // Lower sentinel + _lows[0] = -std::numeric_limits<double>::infinity(); + + // Copy a sorted vector + for (size_t i = 0; i < lows.size(); i++) { + _lows[i+1] = lows[i]; + } + + _max = lows.size() + 1; + + // Upper sentinel + _lows[_max] = std::numeric_limits<double>::infinity(); + } + + // Linear search in the forward direction + size_t _linsearch_forward(double *arr, double key, size_t n) const { + for (size_t i=0; i<n; i++) + if (arr[i] > key) return i; + return n; + } + + // Linear search in the backward direction + size_t _linsearch_backward(double *arr, double key, size_t n) const { + arr -= n; + for (size_t i=0; i < n; i++) + if (arr[n - 1 - i] <= key) return i; + return n; + } + + // Bisection search, adapted from C++ std lib implementation + size_t _bisect(double *arr, + const double key, size_t min, size_t max) const { + size_t len = max - min; + + size_t middle; + while (len > BISECT_LINEAR_THRESHOLD) { + size_t half = len >> 1; + middle = min + half; + if (arr[middle] < key) { + min = middle + 1; + len = len - half - 1; + } else + len = half; + } + return min + _linsearch_forward(arr + min, key, BISECT_LINEAR_THRESHOLD); + } + + protected: + boost::shared_ptr<Estimator> _est; + boost::shared_array<double> _lows; + size_t _max; + }; + } +} + +#endif Added: trunk/include/YODA/Utils/fastlog.h ============================================================================== --- /dev/null 00:00:00 1970 (empty, because file is newly added) +++ trunk/include/YODA/Utils/fastlog.h Tue Aug 7 02:08:15 2012 (r523) @@ -0,0 +1,71 @@ +/*=====================================================================* + * Copyright (C) 2011 Paul Mineiro * + * All rights reserved. * + * * + * Redistribution and use in source and binary forms, with * + * or without modification, are permitted provided that the * + * following conditions are met: * + * * + * * Redistributions of source code must retain the * + * above copyright notice, this list of conditions and * + * the following disclaimer. * + * * + * * Redistributions in binary form must reproduce the * + * above copyright notice, this list of conditions and * + * the following disclaimer in the documentation and/or * + * other materials provided with the distribution. * + * * + * * Neither the name of Paul Mineiro nor the names * + * of other contributors may be used to endorse or promote * + * products derived from this software without specific * + * prior written permission. * + * * + * THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND * + * CONTRIBUTORS "AS IS" AND ANY EXPRESS OR IMPLIED WARRANTIES, * + * INCLUDING, BUT NOT LIMITED TO, THE IMPLIED WARRANTIES * + * OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE * + * ARE DISCLAIMED. IN NO EVENT SHALL THE COPYRIGHT OWNER * + * OR CONTRIBUTORS BE LIABLE FOR ANY DIRECT, INDIRECT, * + * INCIDENTAL, SPECIAL, EXEMPLARY, OR CONSEQUENTIAL DAMAGES * + * (INCLUDING, BUT NOT LIMITED TO, PROCUREMENT OF SUBSTITUTE * + * GOODS OR SERVICES; LOSS OF USE, DATA, OR PROFITS; OR * + * BUSINESS INTERRUPTION) HOWEVER CAUSED AND ON ANY THEORY OF * + * LIABILITY, WHETHER IN CONTRACT, STRICT LIABILITY, OR TORT * + * (INCLUDING NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY * + * OUT OF THE USE OF THIS SOFTWARE, EVEN IF ADVISED OF THE * + * POSSIBILITY OF SUCH DAMAGE. * + * * + * Contact: Paul Mineiro <paul at mineiro.com> * + *=====================================================================*/ + +#ifndef __FAST_LOG_H_ +#define __FAST_LOG_H_ + +#include <stdint.h> + +namespace YODA { + namespace Utils { + static inline float + fastlog2 (float x) + { + union { float f; uint32_t i; } vx = { x }; + union { uint32_t i; float f; } mx = { (vx.i & 0x007FFFFF) | 0x3f000000 }; + float y = vx.i; + y *= 1.1920928955078125e-7f; + + return y - 124.22551499f + - 1.498030302f * mx.f + - 1.72587999f / (0.3520887068f + mx.f); + } + + static inline float + fasterlog2 (float x) + { + union { float f; uint32_t i; } vx = { x }; + float y = vx.i; + y *= 1.1920928955078125e-7f; + return y - 126.94269504f; + } + } +} +#endif // __FAST_LOG_H_ Modified: trunk/pyext/Makefile.am ============================================================================== --- trunk/pyext/Makefile.am Mon Aug 6 14:49:55 2012 (r522) +++ trunk/pyext/Makefile.am Tue Aug 7 02:08:15 2012 (r523) @@ -4,7 +4,6 @@ CYTHONFLAGS = @CYTHONFLAGS@ all-local: - $(PYTHON) setup.py build $(PYTHON) setup.py install --install-lib=build/ install-exec-local: Modified: trunk/pyext/setup.py.in ============================================================================== --- trunk/pyext/setup.py.in Mon Aug 6 14:49:55 2012 (r522) +++ trunk/pyext/setup.py.in Tue Aug 7 02:08:15 2012 (r523) @@ -27,11 +27,6 @@ library_dirs=[srcdir, os.path.join(srcdir, '.libs')], libraries=['stdc++', 'YODA']) -extns = [ext('util'), - ext('core', statics=[statics], - depends=glob('yoda/include/*.pyx') - ), - ] # Make the templates make_templates('Axis1D_BIN1D_DBN', @@ -45,6 +40,15 @@ make_templates('Bin1D_DBN', DBN=('Dbn1D', 'Dbn2D')) make_templates('Bin2D_DBN', DBN=('Dbn2D', 'Dbn3D')) +extns = [ext('util'), + ext('core', statics=[statics], + depends=glob('yoda/include/*.pyx') + ), + ext('axis', statics=[statics], + depends=glob('yoda/include/Axis*.pyx') + ), + ] + setup( name = PKGNAME, ext_modules = extns, Added: trunk/pyext/yoda/axis.pyx ============================================================================== --- /dev/null 00:00:00 1970 (empty, because file is newly added) +++ trunk/pyext/yoda/axis.pyx Tue Aug 7 02:08:15 2012 (r523) @@ -0,0 +1,17 @@ +cimport yoda.declarations as c +cimport yoda.util as util +import yoda.util as util + +from cython.operator cimport dereference as deref, preincrement as preinc +from libcpp.string cimport string +from libcpp.vector cimport vector +from libcpp.pair cimport pair +from libcpp.map cimport map + +# Pure python imports +from itertools import repeat, imap +from operator import attrgetter +from yoda import * + +include "include/Errors.pyx" +include "include/Axis1D.pxi" Modified: trunk/pyext/yoda/core.pyx ============================================================================== --- trunk/pyext/yoda/core.pyx Mon Aug 6 14:49:55 2012 (r522) +++ trunk/pyext/yoda/core.pyx Tue Aug 7 02:08:15 2012 (r523) @@ -20,8 +20,6 @@ include "include/Point3D.pyx" include "include/Bin.pyx" -include "include/Axis1D.pxi" - include "include/Bin1D.pxi" include "include/HistoBin1D.pyx" Modified: trunk/pyext/yoda/declarations.pxd ============================================================================== --- trunk/pyext/yoda/declarations.pxd Mon Aug 6 14:49:55 2012 (r522) +++ trunk/pyext/yoda/declarations.pxd Tue Aug 7 02:08:15 2012 (r523) @@ -767,47 +767,6 @@ operator / (Histo2D) # Histo2D }}} -"""INTERNAL CLASSES -- FOR TESTING PURPOSES""" - -# Axis1D {{{ -cdef extern from "YODA/Axis1D.h" namespace "YODA": - cdef cppclass Axis1D[BIN1D, DBN]: - Axis1D() except+ err - Axis1D(vector[double] binedges) except+ err - Axis1D(size_t, double, double) except+ err - Axis1D(vector[BIN1D] bins) except+ err - size_t numBins() except+ err - vector[BIN1D] &bins() - double lowEdge() except+ err - double highEdge() except+ err - long getBinIndex(double) - void reset() - DBN &totalDbn() - DBN &underflow() - DBN &overflow() -# Axis1D }}} - -# Axis2D {{{ -cdef extern from "YODA/Axis2D.h" namespace "YODA": - cdef cppclass Axis2D[BIN2D, DBN]: - Axis2D() except+ err - Axis2D(vector[double] xedges, vector[double] yedges) except+ err - Axis2D(size_t, pair[double, double], - size_t, pair[double, double]) except+ err - Axis2D(vector[BIN2D] bins) 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() except +err - DBN outflow(int, int) except +err - -# Axis2D }}} - # Streams {{{ cdef extern from "<sstream>" namespace "std": cdef cppclass ostringstream: @@ -840,3 +799,24 @@ cdef extern from "YODA/WriterAIDA.h" namespace "YODA": Writer& WriterAIDA_create "YODA::WriterAIDA::create" () # Streams }}} + +# Axis1D {{{ +cdef extern from "YODA/Axis1D.h" namespace "YODA": + cdef cppclass Axis1D[BIN1D, DBN]: + Axis1D() except+ err + Axis1D(vector[double] binedges) except+ err + Axis1D(size_t, double, double) except+ err + Axis1D(vector[BIN1D] bins) except+ err + void addBin(double, double) except+ err + size_t numBins() except+ err + vector[BIN1D] &bins() + double lowEdge() except+ err + double highEdge() except+ err + long getBinIndex(double) + void reset() + DBN &totalDbn() + DBN &underflow() + DBN &overflow() + void eraseBin(size_t index) except+ err + void mergeBins(size_t, size_t) except+ err +# Axis1D }}} Modified: trunk/pyext/yoda/include/Axis1D_BIN1D_DBN.pyx ============================================================================== --- trunk/pyext/yoda/include/Axis1D_BIN1D_DBN.pyx Mon Aug 6 14:49:55 2012 (r522) +++ trunk/pyext/yoda/include/Axis1D_BIN1D_DBN.pyx Tue Aug 7 02:08:15 2012 (r523) @@ -2,19 +2,22 @@ # it be a user facing class... it's merely there for tests) cdef class Axis1D_${BIN1D}_${DBN}(util.Base): - def __init__(self, size_t nbins, double lower, double upper): + def __init__(self): util.set_owned_ptr( - self, new c.Axis1D[c.${BIN1D}, c.${DBN}]( - nbins, lower, upper)) + self, new c.Axis1D[c.${BIN1D}, c.${DBN}]()) def __len__(self): return self._Axis1D().bins().size() + def __getitem__(self, py_ix): cdef size_t i = util.pythonic_index(py_ix, self._Axis1D().bins().size()) return util.new_borrowed_cls( ${BIN1D}, & self._Axis1D().bins().at(i), self) + def addBin(self, a, b): + self._Axis1D().addBin(a, b) + def __repr__(self): return "<Axis1D>" @@ -36,8 +39,17 @@ def reset(self): self._Axis1D().reset() - def binByCoord(self, x): - return self[self._Axis1D().getBinIndex(x)] + def eraseBin(self, i): + self._Axis1D().eraseBin(i) + + def getBinIndex(self, x): + return self._Axis1D().getBinIndex(x) + + def mergeBins(self, a, b): + self._Axis1D().mergeBins(a, b) + + #def binByCoord(self, x): + # return self[self._Axis1D().getBinIndex(x)] # BOILERPLATE STUFF cdef inline c.Axis1D[c.${BIN1D}, c.${DBN}] *_Axis1D(self) except NULL:
More information about the yoda-svn mailing list |