From 2e5294e1fa0d89fd1c01fb0619060ccb8a63d508 Mon Sep 17 00:00:00 2001 From: Hans Dembinski Date: Sat, 16 Jul 2016 22:07:24 -0400 Subject: [PATCH] super fast regular_axis and code cleanup --- build/CMakeLists.txt | 7 +-- include/boost/histogram/axis.hpp | 91 ++++++++++++++------------- include/boost/histogram/histogram.hpp | 32 +++------- src/axis.cpp | 91 +++------------------------ test/check/speed_cpp.cpp | 79 ++++++++++------------- 5 files changed, 99 insertions(+), 201 deletions(-) diff --git a/build/CMakeLists.txt b/build/CMakeLists.txt index c33ba8d7..cdbc21f7 100644 --- a/build/CMakeLists.txt +++ b/build/CMakeLists.txt @@ -6,7 +6,6 @@ list(APPEND CMAKE_MODULE_PATH ${PROJECT_SOURCE_DIR}) # setup build option(BUILD_CHECKS "Build auxillary checks" OFF) -option(COVERAGE "Build with coverage instrumentation" OFF) if(${CMAKE_BUILD_TYPE}) STRING(TOLOWER ${CMAKE_BUILD_TYPE} CMAKE_BUILD_TYPE) @@ -24,11 +23,11 @@ include_directories(../include ${Boost_INCLUDE_DIRS}) if(CMAKE_BUILD_TYPE STREQUAL "debug") add_definitions(-O0 -g) - set(CMAKE_CXX_FLAGS ${CMAKE_CXX_FLAGS} --coverage) - message(STATUS "Build type: COVERAGE [optimizations off]") + message(STATUS "Build type: DEBUG [optimizations off]") elseif(CMAKE_BUILD_TYPE STREQUAL "cov") add_definitions(-O0 -g) - message(STATUS "Build type: DEBUG [optimizations off]") + set(CMAKE_CXX_FLAGS ${CMAKE_CXX_FLAGS} --coverage) + message(STATUS "Build type: COVERAGE [optimizations off]") elseif(CMAKE_BUILD_TYPE STREQUAL "perf") message(STATUS "Build type: PERF [optimizations on]") add_definitions(-O3 -g -fno-omit-frame-pointer) diff --git a/include/boost/histogram/axis.hpp b/include/boost/histogram/axis.hpp index 55752f01..560d8c5d 100644 --- a/include/boost/histogram/axis.hpp +++ b/include/boost/histogram/axis.hpp @@ -9,7 +9,6 @@ #include #include -#include #include #include #include @@ -31,25 +30,28 @@ namespace histogram { class axis_base { public: ///Returns the number of bins. - inline int bins() const { return size_ > 0 ? size_ : -size_; } + inline int bins() const { return size_ ; } ///Returns whether overflow and underflow bins will be added in the histogram. - inline bool uoflow() const { return size_ > 0; } + inline bool uoflow() const { return uoflow_; } ///Change the label of an axis (not implemented for category_axis). const std::string& label() const { return label_; } ///Returns the axis label, which is a name or description (not implemented for category_axis). void label(const std::string& label) { label_ = label; } protected: - axis_base(int, const std::string&, bool); + axis_base(unsigned, const std::string&, bool); - axis_base() : size_(0) {} - axis_base(const axis_base&); - axis_base& operator=(const axis_base&); + axis_base() = default; + axis_base(const axis_base&) = default; + axis_base(axis_base&&) = default; + axis_base& operator=(const axis_base&) = default; + axis_base& operator=(axis_base&&) = default; bool operator==(const axis_base&) const; private: int size_; + bool uoflow_; std::string label_; template @@ -86,24 +88,27 @@ public: * \param label description of the axis * \param uoflow add underflow and overflow bins to the histogram for this axis or not */ - regular_axis(int n, double min, double max, + regular_axis(unsigned n, double min, double max, const std::string& label = std::string(), bool uoflow = true); - regular_axis() : min_(0), range_(0) {} - regular_axis(const regular_axis&); - regular_axis& operator=(const regular_axis&); + regular_axis() : min_(0), delta_(0) {} + regular_axis(const regular_axis&) = default; + regular_axis(regular_axis&&) = default; + regular_axis& operator=(const regular_axis&) = default; + regular_axis& operator=(regular_axis&&) = default; - ///Returns the bin index for the passed argument. - inline int index(double x) const { - const double z = (x - min_) / range_; - return algorithm::clamp(static_cast(floor(z * bins())), -1, bins()); + ///Returns the bin index for the passed argument (optimized code). + inline int index(double x) const + { + const double z = (x - min_) / delta_; + return z >= 0.0 ? std::min(static_cast(z), bins()) : -1; } double operator[](int idx) const; bool operator==(const regular_axis&) const; private: - double min_, range_; + double min_, delta_; template friend void serialize(Archive& ar, regular_axis & axis ,unsigned version); @@ -113,17 +118,19 @@ private: class polar_axis: public axis_base, public real_axis { public: /** Constructor - \param n number of bins - \param start starting phase of the angle - \param label description of the axis - */ - explicit - polar_axis(int n, double start = 0.0, + * \param n number of bins + * \param start starting phase of the angle + * \param label description of the axis + */ + explicit + polar_axis(unsigned n, double start = 0.0, const std::string& label = std::string()); polar_axis() : start_(0) {} - polar_axis(const polar_axis&); - polar_axis& operator=(const polar_axis&); + polar_axis(const polar_axis&) = default; + polar_axis(polar_axis&&) = default; + polar_axis& operator=(const polar_axis&) = default; + polar_axis& operator=(polar_axis&&) = default; ///Returns the bin index for the passed argument. inline int index(double x) const { @@ -147,9 +154,9 @@ class variable_axis : public axis_base, public real_axis { public: /** Constructor * - * \param x An axis for real-valued data and bins of varying width. Binning is a O(log(N)) operation. If speed matters and the problem domain allows it, prefer a regular_axis. + * \param x An axis for real-valued data and bins of varying width. Binning is a O(log(N)) operation. If speed matters and the problem domain allows it, prefer a regular_axis. * \param label description of the axis - * \param uoflow add underflow and overflow bins to the histogram for this axis or not + * \param uoflow add under-/overflow bins for this axis or not */ explicit variable_axis(const std::initializer_list& x, @@ -187,7 +194,7 @@ public: double operator[](int idx) const; bool operator==(const variable_axis&) const; private: - boost::scoped_array x_; + std::unique_ptr x_; // smaller size compared to std::vector template friend void serialize(Archive& ar, variable_axis & axis, unsigned version); @@ -207,15 +214,16 @@ public: category_axis(const std::initializer_list& categories); category_axis() {} - category_axis(const category_axis&); + category_axis(const category_axis&) = default; category_axis(category_axis&&) = default; - category_axis& operator=(const category_axis&); + category_axis& operator=(const category_axis&) = default; category_axis& operator=(category_axis&&) = default; inline int bins() const { return categories_.size(); } ///Returns the bin index for the passed argument. - inline int index(double x) const { return static_cast(x + 0.5); } + inline int index(double x) const + { return static_cast(x + 0.5); } ///Returns the category for the bin index. const std::string& operator[](int idx) const { return categories_[idx]; } @@ -247,8 +255,10 @@ public: bool uoflow = true); integer_axis() : min_(0) {} - integer_axis(const integer_axis&); - integer_axis& operator=(const integer_axis&); + integer_axis(const integer_axis&) = default; + integer_axis(integer_axis&&) = default; + integer_axis& operator=(const integer_axis&) = default; + integer_axis& operator=(integer_axis&&) = default; ///Returns the bin index for the passed argument. inline int index(double x) const @@ -274,20 +284,13 @@ namespace detail { template void operator()(const A& a) { int j = std::is_same::value ? a.index(in) : in; - const int bins = a.bins(); - const int range = bins + 2 * a.uoflow(); - // the following three lines worout for any overflow setting - j += (j < 0) * (bins + 2); // wrap around if in < 0 - if (j < range) { - out += j * stride; - stride *= range; - } else { - out = std::size_t(-1); - stride = 0; - } + const int shape = a.bins() + 2 * a.uoflow(); + j += (j < 0) * (a.bins() + 2); // wrap around if j < 0 + out += j * stride; +#pragma GCC diagnostic ignored "-Wstrict-overflow" + stride *= (j < shape) * shape; // stride == 0 indicates out-of-range } }; - typedef linearize_t linearize; typedef linearize_t linearize_x; } diff --git a/include/boost/histogram/histogram.hpp b/include/boost/histogram/histogram.hpp index 64ed3d4a..360476de 100644 --- a/include/boost/histogram/histogram.hpp +++ b/include/boost/histogram/histogram.hpp @@ -10,7 +10,10 @@ #include #include #include +#include #include +#include +#include #include #include #include @@ -21,29 +24,10 @@ namespace boost { namespace histogram { - template - struct static_storage_policy { - std::vector data_; - typedef T value_t; - typedef T variance_t; - std::size_t size() const { return data_.size(); } - void allocate(std::size_t n) { data_.resize(n, 0); } - void increase(std::size_t i) { ++data_[i]; } - value_t value(std::size_t i) const { return data_[i]; } - variance_t variance(std::size_t i) const { return data_[i]; } - bool operator==(const static_storage_policy& other) const - { return data_ == other.data_; } - void operator+=(const static_storage_policy& other) - { - for (std::size_t i = 0, n = size(); i < n; ++i) - data_[i] += other.data_[i]; - } - }; - ///Use dynamic dimension constexpr unsigned Dynamic = 0; - template > + template > class histogram_t : private StoragePolicy { public: @@ -95,7 +79,7 @@ namespace histogram { "number of arguments does not match histogram dimension"); detail::linearize_x lin; fill_impl(lin, std::forward(args)...); - if (lin.out < size()) + if (lin.stride) StoragePolicy::increase(lin.out); } @@ -106,7 +90,7 @@ namespace histogram { "number of arguments does not match histogram dimension"); detail::linearize lin; index_impl(lin, std::forward(args)...); - if (lin.out >= size()) + if (lin.stride == 0) throw std::out_of_range("invalid index"); return StoragePolicy::value(lin.out); } @@ -118,7 +102,7 @@ namespace histogram { "number of arguments does not match histogram dimension"); detail::linearize lin; index_impl(lin, std::forward(args)...); - if (lin.out >= size()) + if (lin.stride == 0) throw std::out_of_range("invalid index"); return StoragePolicy::variance(lin.out); } @@ -158,7 +142,7 @@ namespace histogram { } private: - std::array axes_; + std::array axes_; std::size_t field_count() const { diff --git a/src/axis.cpp b/src/axis.cpp index 7160ca1e..b74ad525 100644 --- a/src/axis.cpp +++ b/src/axis.cpp @@ -29,62 +29,27 @@ namespace { } } -axis_base::axis_base(int size, +axis_base::axis_base(unsigned size, const std::string& label, bool uoflow) : size_(size), + uoflow_(uoflow), label_(label) -{ - if (size <= 0) - throw std::logic_error("size must be positive"); - if (!uoflow) size_ = -size_; -} - -axis_base::axis_base(const axis_base& o) : - size_(o.size_), - label_(o.label_) {} -axis_base& -axis_base::operator=(const axis_base& o) -{ - if (this != &o) { - size_ = o.size_; - label_ = o.label_; - } - return *this; -} - bool axis_base::operator==(const axis_base& o) const { return size_ == o.size_ && label_ == o.label_; } -regular_axis::regular_axis(int n, double min, double max, +regular_axis::regular_axis(unsigned n, double min, double max, const std::string& label, bool uoflow): axis_base(n, label, uoflow), min_(min), - range_(max - min) + delta_((max - min) / n) { if (min >= max) throw std::logic_error("regular_axis: min must be less than max"); } -regular_axis::regular_axis(const regular_axis& o) : - axis_base(o), - min_(o.min_), - range_(o.range_) -{} - -regular_axis& -regular_axis::operator=(const regular_axis& o) -{ - if (this != &o) { - axis_base::operator=(o); - min_ = o.min_; - range_ = o.range_; - } - return *this; -} - double regular_axis::operator[](int idx) const @@ -94,7 +59,7 @@ regular_axis::operator[](int idx) if (idx > bins()) return std::numeric_limits::infinity(); const double z = double(idx) / bins(); - return (1.0 - z) * min_ + z * (min_ + range_); + return (1.0 - z) * min_ + z * (min_ + delta_ * bins()); } bool @@ -102,30 +67,15 @@ regular_axis::operator==(const regular_axis& o) const { return axis_base::operator==(o) && min_ == o.min_ && - range_ == o.range_; + delta_ == o.delta_; } -polar_axis::polar_axis(int n, double start, +polar_axis::polar_axis(unsigned n, double start, const std::string& label) : axis_base(n, label, false), start_(start) {} -polar_axis::polar_axis(const polar_axis& o) : - axis_base(o), - start_(o.start_) -{} - -polar_axis& -polar_axis::operator=(const polar_axis& o) -{ - if (this != &o) { - axis_base::operator=(o); - start_ = o.start_; - } - return *this; -} - double polar_axis::operator[](int idx) const @@ -186,18 +136,6 @@ category_axis::category_axis(const std::initializer_list& c) : categories_(c) {} -category_axis::category_axis(const category_axis& o) : - categories_(o.categories_) -{} - -category_axis& -category_axis::operator=(const category_axis& o) -{ - if (this != &o) - categories_ = o.categories_; - return *this; -} - bool category_axis::operator==(const category_axis& o) const { return categories_ == o.categories_; } @@ -209,21 +147,6 @@ integer_axis::integer_axis(int min, int max, min_(min) {} -integer_axis::integer_axis(const integer_axis& a) : - axis_base(a), - min_(a.min_) -{} - -integer_axis& -integer_axis::operator=(const integer_axis& o) -{ - if (this != &o) { - axis_base::operator=(o); - min_ = o.min_; - } - return *this; -} - bool integer_axis::operator==(const integer_axis& o) const { diff --git a/test/check/speed_cpp.cpp b/test/check/speed_cpp.cpp index 6edb2fe4..36c39618 100644 --- a/test/check/speed_cpp.cpp +++ b/test/check/speed_cpp.cpp @@ -6,54 +6,44 @@ #include #include -#include -#include +#include #include #include -#include +#include #include #include -using namespace std; using namespace boost::histogram; -template -struct rng { - boost::random::mt19937 r; - D d; - rng(double a, double b) : d(a, b) {} - double operator()() { return d(r); } -}; - -vector random_array(unsigned n, int type) { - using namespace boost::random; - std::vector result; - switch (type) { - case 0: - std::generate_n(std::back_inserter(result), n, rng >(0.0, 1.0)); - break; - case 1: - std::generate_n(std::back_inserter(result), n, rng >(0.0, 0.3)); - break; +std::vector random_array(unsigned n, int type) { + std::vector result(n); + std::default_random_engine gen(1); + if (type) { // type == 1 + std::normal_distribution<> d(0.0, 0.3); + for (auto& x : result) + x = d(gen); + } + else { // type == 0 + std::uniform_real_distribution<> d(0.0, 1.0); + for (auto& x: result) + x = d(gen); } return result; } void compare_1d(unsigned n, int distrib) { - vector r = random_array(n, distrib); + auto r = random_array(n, distrib); - double best_boost = std::numeric_limits::max(); + auto best_boost = std::numeric_limits::max(); for (unsigned k = 0; k < 10; ++k) { - histogram h(regular_axis(100, 0, 1)); - t = clock(); + auto h = histogram(regular_axis(100, 0, 1)); + auto t = clock(); for (unsigned i = 0; i < n; ++i) h.fill(r[i]); t = clock() - t; best_boost = std::min(best_boost, double(t) / CLOCKS_PER_SEC); - // printf("root %g this %g\n", hroot.GetSum(), h.sum()); - assert(hroot.GetSum() == h.sum()); } printf("1D\n"); @@ -62,19 +52,18 @@ void compare_1d(unsigned n, int distrib) void compare_3d(unsigned n, int distrib) { - vector r = random_array(3 * n, distrib); + auto r = random_array(3 * n, distrib); - double best_boost = std::numeric_limits::max(); + auto best_boost = std::numeric_limits::max(); for (unsigned k = 0; k < 10; ++k) { - histogram h(regular_axis(100, 0, 1), - regular_axis(100, 0, 1), - regular_axis(100, 0, 1)); - t = clock(); + auto h = histogram(regular_axis(100, 0, 1), + regular_axis(100, 0, 1), + regular_axis(100, 0, 1)); + auto t = clock(); for (unsigned i = 0; i < n; ++i) h.fill(r[3 * i], r[3 * i + 1], r[3 * i + 2]); t = clock() - t; best_boost = std::min(best_boost, double(t) / CLOCKS_PER_SEC); - assert(hroot.GetSum() == h.sum()); } printf("3D\n"); @@ -83,24 +72,24 @@ void compare_3d(unsigned n, int distrib) void compare_6d(unsigned n, int distrib) { - vector r = random_array(6 * n, distrib); + auto r = random_array(6 * n, distrib); - double best_boost = std::numeric_limits::max(); + auto best_boost = std::numeric_limits::max(); for (unsigned k = 0; k < 10; ++k) { double x[6]; - histogram h(regular_axis(10, 0, 1), - regular_axis(10, 0, 1), - regular_axis(10, 0, 1), - regular_axis(10, 0, 1), - regular_axis(10, 0, 1), - regular_axis(10, 0, 1)); + auto h = histogram(regular_axis(10, 0, 1), + regular_axis(10, 0, 1), + regular_axis(10, 0, 1), + regular_axis(10, 0, 1), + regular_axis(10, 0, 1), + regular_axis(10, 0, 1)); - t = clock(); + auto t = clock(); for (unsigned i = 0; i < n; ++i) { for (unsigned k = 0; k < 6; ++k) x[k] = r[6 * i + k]; - h.fill(x, x+6); + h.fill(x[0], x[1], x[2], x[3], x[4], x[5]); } t = clock() - t; best_boost = std::min(best_boost, double(t) / CLOCKS_PER_SEC);