From a40b884a2bf12e3d688eb05b6273b890f0f1dc97 Mon Sep 17 00:00:00 2001 From: Ryan Date: Thu, 8 Jun 2023 21:08:47 -0400 Subject: [PATCH 01/10] getting basic idea in --- include/boost/histogram/axis/piecewise.hpp | 639 +++++++++++++++++++++ 1 file changed, 639 insertions(+) create mode 100644 include/boost/histogram/axis/piecewise.hpp diff --git a/include/boost/histogram/axis/piecewise.hpp b/include/boost/histogram/axis/piecewise.hpp new file mode 100644 index 00000000..7c4fdabb --- /dev/null +++ b/include/boost/histogram/axis/piecewise.hpp @@ -0,0 +1,639 @@ +// Copyright 2015-2018 Hans Dembinski +// +// Distributed under the Boost Software License, Version 1.0. +// (See accompanying file LICENSE_1_0.txt +// or copy at http://www.boost.org/LICENSE_1_0.txt) + +#ifndef BOOST_HISTOGRAM_AXIS_PIECEWISE_HPP +#define BOOST_HISTOGRAM_AXIS_PIECEWISE_HPP + +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include + +namespace boost { +namespace histogram { + +namespace axis { + +/**Pieces*** + + Each piece of a piecewise axis has a start point x0 and an end point x1. + + |--------------------------------------| + x0 x1 + + Each piece has n bins. There are three options for varying bin sizes. + + Constant |------|------|------|------| + (bi+1 = bi) bi bi+1 + + Multiply |---|-------|---------------| + (bi+1 = bi * constant) bi bi+1 + + Add |---|-----|-------|---------| + (bi+1 = bi + constant) bi bi+1 + +***Creating a Piece*** + + To create a piece, one starts at a point x0. Starting with a bin size b0, bins are + added according to the bin spacing option until a termination condition is met. The + termination condition can be either a bin number n_bins or stopping point x_stop. + These two possibilities are shown below. + + Add n bins + Start at x0 = 1. Start with a bin size of b0 = 1. Double the bin size until 4 bins + are added. + |-|---|-------|---------------| + 0 2 4 6 8 10 12 14 16 + + Add until the piece contains x_stop + Start at x0 = 1. Start with a bin size of b0 = 3. Add bins of constant size until + the piece contains x_stop = 14. + |-----|-----|-----|-----|--*--| + 0 2 4 6 8 10 12 14 16 + +***Extrapolation and Attachment*** + + Pieces can be added to the right or left side of an existing piece. Pieces are + added with one of two approaches: extrapolation or attachment. Examples of each are + shown below. + Extrapolate 2 bins to the right side, doubling the bin size each time. + 0 2 4 6 8 + |---|---|---|---| + <- - - New piece - - -> + |---|---|---|---|-------|---------------| + 0 2 4 6 8 10 12 14 16 18 20 + Attach a new piece with a uniform bin spacing to the right side. + 0 2 4 6 8 + |---|---|---|---| + <- - - New piece - - -> + |---|---|---|---|-|-|-|-|-|-|-|-|-|-|-|-| + 0 2 4 6 8 10 12 14 16 18 20 + +***Syntax*** + +using P = piecewise; + +auto pa = P( // Create a new piece + P::x0(1), // Start at x0 = 1 + P::b0(1), // Start with a bin size of b0 = 1 + P::bin_trans::uniform(1), // Add bins of constant size 1 + P::stop::n_bins(4)); // Add 4 bins + +pa.add_right( // Add to the right side + P::bin_trans::multiply(2), // Extrapolate, doubling the bin size each time + P::stop::x(16)); // Stop when the piece contains x = 16 + +pa.add_left( // Add to the left side + P::b0(1), // Start with a bin size of b0 = 1 + P::bin_trans::add(1), // Add bins of constant size 1 + P::stop::x(4)); // Stop when the piece contains x = 4 +*/ + +/// Struct for starting point +class x0 { +public: + x0(double value) : value_{value} { + if (!std::isfinite(value)) + BOOST_THROW_EXCEPTION(std::invalid_argument("x0 must be finite")); + } + + double value() const { return value_; } + +private: + double value_; +}; + +/// Struct for bin size b0 +class b0 { +public: + b0(double value) : value_{value} { + if (!(0 < value)) BOOST_THROW_EXCEPTION(std::invalid_argument("b0 must be > 0")); + } + + double value() const { return value_; } + +private: + double value_; +}; + +/// Class for bin spacing option and value +class bin_transion { +public: + enum class option { uniform, multiply, add }; + + double calc_next_bin_size(double b) const { + assert(0 < b); + switch (option_) { + case option::uniform: assert(b == value_); return b; + case option::multiply: return b * value_; + case option::add: return b + value_; + } + } + + bin_transion uniform(double value) = bin_transion{option::uniform, value}; + bin_transion multiply(double value) = bin_transion{option::multiply, value}; + bin_transion add(double value) = bin_transion{option::add, value}; + + bool is_uniform() const { return option_ == option::uniform; } + bool is_multiply() const { return option_ == option::multiply; } + bool is_add() const { return option_ == option::add; } + + double value() const { return value_; } + +private: + bin_transion(option option, double value) : option_{option}, value_{value} {} + + option option_; + double value_; +}; + +/// Class for stopping condition and value +class bin_stop { +public: + enum class option { n_bins, x }; + + bin_stop n_bins(double value) = bin_stop{option::n_bins, value}; + bin_stop x(double value) = bin_stop{option::x, value}; + + bool is_n_bins() const { return option_ == option::n_bins; } + bool is_x() const { return option_ == option::x; } + + double value() const { return value_; } + +private: + bin_stop(option option, double value) : option_{option}, value_{value} {} + + option option_; + double value_; +}; + +/// +template +class piece { + using value_type = Value; + using internal_value_type = detail::get_scale_type; + +public: + piece(x0 x0, b0 b0, bin_trans bin_trans, bin_stop stop, bool is_right) + : x0_{x0.value}, b0_{b0.value}, bin_trans_{bin_trans} { + // Assert is handled case + if (stop.is_n_bins()) { + if (stop.n_bins.value <= 0) { + BOOST_THROW_EXCEPTION(std::invalid_argument("n_bins must be > 0")); + } + } else if (stop.is_x()) { + if (stop.x.value <= x0_) { + BOOST_THROW_EXCEPTION(std::invalid_argument("x must be > x0")); + } + } + + // Initialize + initialize(stop); + } + + /// Add bins left + void add_left_bins(index_type n) noexcept { + assert(0 < n); // TODO: remove + n = std::abs(n); + + // Update x0 and bin size + x0_ = reverse(-n); + b0_ = reverse(-n + 1) - x0_; + size_++; + } + + /// Add bins right + void add_bins_right(index_type n) noexcept { + assert(0 < n); // TODO: remove + n = std::abs(n); + + // Update x0 and bin size + xN_ = reverse(size() + n); + b0_ = xN_ - reverse(size() + n - 1); + size_++; + } + + /// Shifts the axis + void shift_axis(index_type n) { + if (n < 0) { + // Shift axis to the left by n bins + const auto stop = min_ + delta_; + min_ += n * (delta_ / size()); + delta_ = stop - min_; + size_ -= n; + + } else if (0 < n) { + // Shift axis to the right by n bins + delta_ /= size(); + delta_ *= size() + n; + size_ += n; + + x0_ += n * b0_; + } + } + + template + T forward(T x) const noexcept { + // Runs in hot loop, please measure impact of changes + switch (option_) { + case option::uniform: return forward_uniform(x); + case option::multiply: return forward_multiply(x); + case option::add: return forward_add(x); + } + } + + index_type size() const noexcept { return size_; } + +private: + // Initialize + // Have: (x0, b0, n_bins) + // Calc: xN + + // Initialize + // Have: (x0, xN, n_bins) + // Calc: b0 + + // Add right + // Have: (x0, b0, n_bins++) + // Update: xN + + // Add left + // Have: (xN, b0, n_bins++) + // Calc: b0 + // Update: x0 + + void initialize(bin_stop stop) { + xN_ = x0_; + size_ = 0; + double b = b0_; + + while (!see_if_need_more_bins()) { + b = bin_trans_.calc_next_bin_size(b); + xN_ += b; + ++size_; + } + } + + bool see_if_need_more_bins() { + if (stop_.is_n_bins()) { + return size_ >= stop_.n_bins.value; + } else if (stop_.is_x()) { + return xN_ >= stop_.x.value; + } + } + + // Uniform spacing + // + // b0 b0 b0 b0 + // N=0 N=1 N=2 N=3 N=4 + // |------|------|------|------| + // x0 x4 + // + // The sequence + // x0 = x0 + // x1 = x0 + b0 + // x2 = x0 + b0 * 2 + // x3 = x0 + b0 * 3 + // ... + // x = x0 + b0 * N + // + // Formula for N in terms of x + // N = (x - x0) / b0 + // + template + interval_value_type forward_linear(T x) const noexcept { + // Runs in hot loop, please measure impact of changes + return ((x / unit_type{}) - x0_) / b0_; + } + + // Multiply spacing + // + // b0 b0*r b0*r² b0*r³ + // N=0 N=1 N=2 N=3 N=4 + // |---|-----|-----------|------------------| + // x0 x4 + // + // The sequence (for some value of ψ) + // x0 = 1 * b0 / (r - 1) + ψ + // x1 = r * b0 / (r - 1) + ψ + // x2 = r² * b0 / (r - 1) + ψ + // x3 = r³ * b0 / (r - 1) + ψ + // ... + // x = rᴺ * b0 / (r - 1) + ψ + // + // Note: the first bin spacing is b0 + // x1 - x0 = r * b0 / (r - 1) - 1 * b0 / (r - 1) + // = (r - 1) * b0 / (r - 1) + // = b0 + // + // Find a formula for N + // x - x0 = b0 / (r - 1) * rᴺ - b0 / (r - 1) * 1 + // = b0 / (r - 1) * (rᴺ - 1) + // + // N in terms of x + // N = log2(((x - x0) / b0) * (r - 1) + 1) / log2(r) + // + template + interval_value_type forward_multiply(T x) const noexcept { + // Runs in hot loop, please measure impact of changes + const auto z = forward_linear(x); // z = (x - x0) / b0 + const auto r = bin_trans_.value(); + return std::log2(z * (r - 1) + 1) / std::log2(r); + } + + // Add spacing + // + // b0 b0+r b0+2r b0+3r + // N=0 N=1 N=2 N=3 N=4 + // |---|-----|-----------|------------------| + // x0 x4 + // + // The sequence + // x0 = x0 + // x1 = x0 + 1 * b0 + // x2 = x0 + 2 * b0 + r + // x3 = x0 + 3 * b0 + r * (1 + 2) + // x4 = x0 + 4 * b0 + r * (1 + 2 + 3) + // ... + // x = x0 + N * b0 + r * (1 + 2 + ... + N-1) + // = x0 + N * b0 + r * (N * (N - 1) / 2) + // + // N in terms of x + // 2 * (x - x0) / b0 = N * (N - 1) + // N² - N - 2 * (x - x0) / b0 = 0 + // Solve quadratic equation and take positive root + // + template + interval_value_type forward_multiply(T x) const noexcept { + // Runs in hot loop, please measure impact of changes + const auto z = forward_linear(x); // z = (x - x0) / b0 + const auto r = bin_trans_.value(); + return (1 + std::sqrt(1 + 8 * z / r)) / 2; // Check copilot suggestion + } + + internal_value_type x0_ { 0 } + double b0_{1}; + bin_trans bin_trans_{bin_trans::uniform{1}}; + + internal_value_type xN_{0} index_type size_{0}; +}; + +/// +class piecewise { +public: + /// Construct a piecewise axis with an initial piece. + piecewise(x0 x0, b0 b0, bin_trans bin_trans, bin_stop stop) { + // + const auto p = piece(x0, b0, bin_trans, stop, true); + v_pieces_.push_back(p); + } + + /// + template + T forward(T x) const noexcept { + // Runs in hot loop, please measure impact of changes + index_type offset = 0; + const int n_pieces = v_pieces_.size(); + for (int j = 0; j < n_pieces; ++j) { + const auto& p = v_pieces_[j]; + if (x < p.x1 || j == n_pieces - 1) return offset + p.index(x); + } + } + + /// + template + index_type index(ValueType x) const noexcept { + // Runs in hot loop, please measure impact of changes + const auto i = forward(x); + if (i < size_) { + if (0 <= i) + return static_cast(i); // 0 <= i < size + else + return -1; // i < 0 + } + // upper edge of last bin is inclusive if overflow bin is not present + if (!options_type::test(option::overflow) && z == 1) return size() - 1; + return size(); // also returned if x is NaN + } + + /// Shifts the axis + void shift_axis(index_type n) { + if (n < 0) { + v_pieces_.front().shift_axis(n); + } else if (0 < n) { + // Get last element of v_pieces + v_pieces_.back().shift_axis(n); + } + } + + /// Returns index and shift (if axis has grown) for the passed argument. + std::pair update(value_type x) noexcept { + assert(options_type::test(option::growth)); + const auto i = forward(x); + if (i < size()) { + if (0 <= i) { + const auto i_int = static_cast(i); + return {i_int, 0}; + } else + (i != -std::numeric_limits::infinity()) { + const auto i_int = static_cast(std::floor(i)); + shift_axis(i_int); + return {0, -i_int}; + } + else { // i is -infinity + return {-1, 0}; + } + } + // i either beyond range, infinite, or NaN + if (z < std::numeric_limits::infinity()) { + const auto i_int = static_cast(i); + const auto n = i_int - size() + 1; + shift_axis(n); + return {i, -n}; + } + // z either infinite or NaN + return {size(), 0}; + } + + add_left(bin_trans bin_trans, bin_stop stop) { + double b0 = get_left_b0(); + add_left(b0, bin_trans, stop); + } + add_right(bin_trans bin_trans, bin_stop stop) { + double b0 = get_right_b0(); + add_right(b0, bin_trans, stop); + } + add_left(b0 b0, bin_trans bin_trans, bin_stop stop) { + double x0 = get_left_x(); + const auto p = piece(x0, b0, bin_trans, stop, false); + v_pieces_.insert(v_pieces_.begin(), p); + } + add_right(b0 b0, bin_trans bin_trans, bin_stop stop) { + double x0 = get_right_x(); + const auto p = piece(x0, b0, bin_trans, stop, true); + v_pieces_.push_back(p); + } + +private: + double get_left_x() const { + assert(!v_pieces_.empty()); + return v_pieces_.front().x0; + } + + double get_right_x() const { + assert(!v_pieces_.empty()); + return v_pieces_.back().x1; + } + + /// Forward transform of external value x. + template + auto forward(T x) const { + return std::pow(x, power); + } + + /// Inverse transform of internal value x. + template + auto inverse(T x) const { + return std::pow(x, 1.0 / power); + } + + std::vector v_pieces_{}; +}; + +/** Axis for unit intervals on the real line. + + The most common binning strategy. Very fast. Binning is a O(1) operation. + + If the axis has an overflow bin (the default), a value on the upper edge of the last + bin is put in the overflow bin. The axis range represents a semi-open interval. + + If the overflow bin is deactivated, then a value on the upper edge of the last bin is + still counted towards the last bin. The axis range represents a closed interval. + + The options `growth` and `circular` are mutually exclusive. + + @tparam Value input value type, must be floating point. + @tparam Transform builtin or user-defined transform type. + @tparam MetaData type to store meta data. + @tparam Options see boost::histogram::axis::option. + */ +template +class unit_regular + : public iterator_mixin>, + protected detail::replace_default, + public metadata_base_t { + // these must be private, so that they are not automatically inherited + using value_type = Value; + using metadata_base = metadata_base_t; + using metadata_type = typename metadata_base::metadata_type; + using options_type = + detail::replace_default; + + using unit_type = detail::get_unit_type; + using internal_value_type = detail::get_scale_type; + +public: + constexpr unit_regular() = default; + + unit_regular(piecewise pw, metadata_type meta = {}, options_type options = {}) + : metadata_base(std::move(meta)), pw_(std::move(pw)) { + // static_asserts were moved here from class scope to satisfy deduction in gcc>=11 + static_assert(std::is_nothrow_move_constructible::value, + "piecewise must be no-throw move constructible"); + static_assert(std::is_nothrow_move_assignable::value, + "piecewise must be no-throw move assignable"); + static_assert(std::is_floating_point::value, + "unit_regular axis requires floating point type"); + static_assert(!(options.test(option::circular) && options.test(option::growth)), + "circular and growth options are mutually exclusive"); + if (size() <= 0) BOOST_THROW_EXCEPTION(std::invalid_argument("bins > 0 required")); + if (!std::isfinite(min_) || !std::isfinite(delta_)) + BOOST_THROW_EXCEPTION( + std::invalid_argument("forward transform of start or stop invalid")); + if (delta_ == 0) + BOOST_THROW_EXCEPTION(std::invalid_argument("range of axis is zero")); + } + + /// Constructor used by algorithm::reduce to shrink and rebin (not for users). + unit_regular(const unit_regular& src, index_type begin, index_type end, unsigned merge) + : unit_regular(src.transform(), (end - begin) / merge, src.value(begin), + src.value(end), src.metadata()) { + assert(false); + assert((end - begin) % merge == 0); + if (options_type::test(option::circular) && !(begin == 0 && end == src.size())) + BOOST_THROW_EXCEPTION(std::invalid_argument("cannot shrink circular axis")); + } + + /// Return index for value argument. + index_type index(value_type x) const noexcept { return piecewise_.index(x); } + + /// Returns index and shift (if axis has grown) for the passed argument. + std::pair update(value_type x) noexcept { + assert(options_type::test(option::growth)); + return piecewise_.update(x); + } + + /// Return value for fractional index argument. + value_type value(real_index_type i) const noexcept { + assert(false); + return static_cast(i); + } + + /// Return bin for index argument. + decltype(auto) bin(index_type idx) const noexcept { + return interval_view(*this, idx); + } + + /// Returns the number of bins, without over- or underflow. + index_type size() const noexcept { return size_; } + + /// Returns the options. + static constexpr unsigned options() noexcept { return options_type::value; } + + // template + // bool operator==(const unit_regular& o) const noexcept { + // return detail::relaxed_equal{}(transform(), o.transform()) && size() == o.size() && + // min_ == o.min_ && delta_ == o.delta_ && + // detail::relaxed_equal{}(this->metadata(), o.metadata()); + // } + // template + // bool operator!=(const unit_regular& o) const noexcept { + // return !operator==(o); + // } + + // template + // void serialize(Archive& ar, unsigned /* version */) { + // ar& make_nvp("transform", static_cast(*this)); + // ar& make_nvp("size", size_); + // ar& make_nvp("meta", this->metadata()); + // ar& make_nvp("min", min_); + // ar& make_nvp("delta", delta_); + // } + +private: + index_type size_{0}; + internal_value_type min_{0} piecewise piecewise_{}; +}; + +} // namespace axis +} // namespace histogram +} // namespace boost + +#endif From 7ff9f6f30401f1eb9c7008fe6f292fe8fa254123 Mon Sep 17 00:00:00 2001 From: Ryan Date: Fri, 9 Jun 2023 19:59:55 -0400 Subject: [PATCH 02/10] created int_resolver and int_resolver_circle --- include/boost/histogram/axis.hpp | 1 + include/boost/histogram/axis/piecewise.hpp | 760 ++++++++++++--------- include/boost/histogram/fwd.hpp | 6 + test/CMakeLists.txt | 1 + test/Jamfile | 1 + test/axis_piecewise_test.cpp | 392 +++++++++++ 6 files changed, 853 insertions(+), 308 deletions(-) create mode 100644 test/axis_piecewise_test.cpp diff --git a/include/boost/histogram/axis.hpp b/include/boost/histogram/axis.hpp index ae42582b..078be5cf 100644 --- a/include/boost/histogram/axis.hpp +++ b/include/boost/histogram/axis.hpp @@ -20,6 +20,7 @@ #include #include #include +#include #include #include #include diff --git a/include/boost/histogram/axis/piecewise.hpp b/include/boost/histogram/axis/piecewise.hpp index 7c4fdabb..aedd12ee 100644 --- a/include/boost/histogram/axis/piecewise.hpp +++ b/include/boost/histogram/axis/piecewise.hpp @@ -26,92 +26,19 @@ #include #include +#include + namespace boost { namespace histogram { namespace axis { -/**Pieces*** - - Each piece of a piecewise axis has a start point x0 and an end point x1. - - |--------------------------------------| - x0 x1 - - Each piece has n bins. There are three options for varying bin sizes. - - Constant |------|------|------|------| - (bi+1 = bi) bi bi+1 - - Multiply |---|-------|---------------| - (bi+1 = bi * constant) bi bi+1 - - Add |---|-----|-------|---------| - (bi+1 = bi + constant) bi bi+1 - -***Creating a Piece*** - - To create a piece, one starts at a point x0. Starting with a bin size b0, bins are - added according to the bin spacing option until a termination condition is met. The - termination condition can be either a bin number n_bins or stopping point x_stop. - These two possibilities are shown below. - - Add n bins - Start at x0 = 1. Start with a bin size of b0 = 1. Double the bin size until 4 bins - are added. - |-|---|-------|---------------| - 0 2 4 6 8 10 12 14 16 - - Add until the piece contains x_stop - Start at x0 = 1. Start with a bin size of b0 = 3. Add bins of constant size until - the piece contains x_stop = 14. - |-----|-----|-----|-----|--*--| - 0 2 4 6 8 10 12 14 16 - -***Extrapolation and Attachment*** - - Pieces can be added to the right or left side of an existing piece. Pieces are - added with one of two approaches: extrapolation or attachment. Examples of each are - shown below. - Extrapolate 2 bins to the right side, doubling the bin size each time. - 0 2 4 6 8 - |---|---|---|---| - <- - - New piece - - -> - |---|---|---|---|-------|---------------| - 0 2 4 6 8 10 12 14 16 18 20 - Attach a new piece with a uniform bin spacing to the right side. - 0 2 4 6 8 - |---|---|---|---| - <- - - New piece - - -> - |---|---|---|---|-|-|-|-|-|-|-|-|-|-|-|-| - 0 2 4 6 8 10 12 14 16 18 20 - -***Syntax*** - -using P = piecewise; - -auto pa = P( // Create a new piece - P::x0(1), // Start at x0 = 1 - P::b0(1), // Start with a bin size of b0 = 1 - P::bin_trans::uniform(1), // Add bins of constant size 1 - P::stop::n_bins(4)); // Add 4 bins - -pa.add_right( // Add to the right side - P::bin_trans::multiply(2), // Extrapolate, doubling the bin size each time - P::stop::x(16)); // Stop when the piece contains x = 16 - -pa.add_left( // Add to the left side - P::b0(1), // Start with a bin size of b0 = 1 - P::bin_trans::add(1), // Add bins of constant size 1 - P::stop::x(4)); // Stop when the piece contains x = 4 -*/ - -/// Struct for starting point -class x0 { +/// Class for starting point +class x_start { public: - x0(double value) : value_{value} { + x_start(double value) : value_{value} { if (!std::isfinite(value)) - BOOST_THROW_EXCEPTION(std::invalid_argument("x0 must be finite")); + BOOST_THROW_EXCEPTION(std::invalid_argument("x_start must be finite")); } double value() const { return value_; } @@ -120,11 +47,12 @@ class x0 { double value_; }; -/// Struct for bin size b0 -class b0 { +/// Class for stopping point +class x_stop { public: - b0(double value) : value_{value} { - if (!(0 < value)) BOOST_THROW_EXCEPTION(std::invalid_argument("b0 must be > 0")); + x_stop(double value) : value_{value} { + if (!std::isfinite(value)) + BOOST_THROW_EXCEPTION(std::invalid_argument("x_stop must be finite")); } double value() const { return value_; } @@ -133,172 +61,240 @@ class b0 { double value_; }; -/// Class for bin spacing option and value -class bin_transion { +/// Class for width_start +class width_start { public: - enum class option { uniform, multiply, add }; - - double calc_next_bin_size(double b) const { - assert(0 < b); - switch (option_) { - case option::uniform: assert(b == value_); return b; - case option::multiply: return b * value_; - case option::add: return b + value_; - } + width_start(double value) : value_{value} { + if (!(0 < value)) + BOOST_THROW_EXCEPTION(std::invalid_argument("width_start must be > 0")); } - bin_transion uniform(double value) = bin_transion{option::uniform, value}; - bin_transion multiply(double value) = bin_transion{option::multiply, value}; - bin_transion add(double value) = bin_transion{option::add, value}; - - bool is_uniform() const { return option_ == option::uniform; } - bool is_multiply() const { return option_ == option::multiply; } - bool is_add() const { return option_ == option::add; } - double value() const { return value_; } private: - bin_transion(option option, double value) : option_{option}, value_{value} {} - - option option_; double value_; }; -/// Class for stopping condition and value -class bin_stop { +/** Class that stores the bin spacing option and constant value + + There are three options for varying bin sizes. + + Uniform |------|------|------|------| + (bi+1 = bi = constant) bi bi+1 + + Multiply |---|-------|---------------| + (bi+1 = bi * constant) bi bi+1 + + Add |---|-----|-------|---------| + (bi+1 = bi + constant) bi bi+1 +*/ +class bin_trans { public: - enum class option { n_bins, x }; + enum class enum_space { e_uniform, e_multiply, e_add }; - bin_stop n_bins(double value) = bin_stop{option::n_bins, value}; - bin_stop x(double value) = bin_stop{option::x, value}; + static bin_trans uniform(double value) { + return bin_trans(enum_space::e_uniform, value); + } + static bin_trans multiply(double value) { + return bin_trans(enum_space::e_multiply, value); + } + static bin_trans add(double value) { return bin_trans(enum_space::e_add, value); } + + bin_trans mirror() const noexcept { + if (is_uniform()) return uniform(value_); + if (is_multiply()) return multiply(1 / value_); + if (is_add()) return add(-value_); + assert(false); // TODO: make unnecessary + } - bool is_n_bins() const { return option_ == option::n_bins; } - bool is_x() const { return option_ == option::x; } + bool is_uniform() const { return spacing_ == enum_space::e_uniform; } + bool is_multiply() const { return spacing_ == enum_space::e_multiply; } + bool is_add() const { return spacing_ == enum_space::e_add; } + enum_space spacing() const { return spacing_; } double value() const { return value_; } private: - bin_stop(option option, double value) : option_{option}, value_{value} {} + bin_trans(enum_space spacing, double value) : spacing_{spacing}, value_{value} {} - option option_; + enum_space spacing_{enum_space::e_uniform}; double value_; }; -/// +/** Class for a piece of a piecewise axis + + Each piece of a piecewise axis has a start point x0 and an end point xN. + It is required that x0 < xN. + + |--------------------------------------| + x0 xN + + Each piece is divided into n bins. The bin spacing is determined by the bin spacing + option and value. + + There are two ways to create a piece: rollout and force. With the rollout method, bins + are added one at a time. If the bin number is positive (n > 0), n bins are added to the + right. If the bin number is negative (n < 0), |n| bins are added to the left. With the + force method, one gives the interval x0 to xN and the number of bins n. The problem + finds the bin width. + +*/ template class piece { using value_type = Value; + using unit_type = detail::get_unit_type; using internal_value_type = detail::get_scale_type; public: - piece(x0 x0, b0 b0, bin_trans bin_trans, bin_stop stop, bool is_right) - : x0_{x0.value}, b0_{b0.value}, bin_trans_{bin_trans} { - // Assert is handled case - if (stop.is_n_bins()) { - if (stop.n_bins.value <= 0) { - BOOST_THROW_EXCEPTION(std::invalid_argument("n_bins must be > 0")); - } - } else if (stop.is_x()) { - if (stop.x.value <= x0_) { - BOOST_THROW_EXCEPTION(std::invalid_argument("x must be > x0")); - } - } + /** Creates a piece with the rollout method - // Initialize - initialize(stop); - } + @param n number of bins (right if positive, left if negative) + @param x_start starting point + @param width_start starting bin width + @param bt bin spacing option and value + + The constructor throws `std::invalid_argument` if the bin number is zero or if the + bin spacing option is invalid. - /// Add bins left - void add_left_bins(index_type n) noexcept { - assert(0 < n); // TODO: remove - n = std::abs(n); + Example: Start at x = 1 with an initial bin size of 1. Rollout 4 bins to the right, + doubling the bin size each time. + |-|---|-------|---------------| + 0 2 4 6 8 10 12 14 16 - // Update x0 and bin size - x0_ = reverse(-n); - b0_ = reverse(-n + 1) - x0_; - size_++; + */ + static piece rollout(int n, x_start x, width_start b0, bin_trans bt) { + assert(n != 0); + if (0 < n) { + return rollout_right(n, x, b0, bt); + } else { + return rollout_left(n, x, b0, bt); + } } - /// Add bins right - void add_bins_right(index_type n) noexcept { - assert(0 < n); // TODO: remove - n = std::abs(n); + /** Creates a piece with the force method - // Update x0 and bin size - xN_ = reverse(size() + n); - b0_ = xN_ - reverse(size() + n - 1); - size_++; - } + @param n number of bins + @param x_start starting point + @param x_stop stopping point + @param bt bin spacing option and value - /// Shifts the axis - void shift_axis(index_type n) { - if (n < 0) { - // Shift axis to the left by n bins - const auto stop = min_ + delta_; - min_ += n * (delta_ / size()); - delta_ = stop - min_; - size_ -= n; + The constructor throws `std::invalid_argument` if: + 1. the bin number is zero, + 2. x_stop <= x_start, or + 3. the bin spacing option is invalid. - } else if (0 < n) { - // Shift axis to the right by n bins - delta_ /= size(); - delta_ *= size() + n; - size_ += n; + Example: Input interval [1, 16] - x0_ += n * b0_; - } + |-----------------------------| + 0 2 4 6 8 10 12 14 16 + + Place 4 bins. Each bin is double the size of the previous bin. + + |-|---|-------|---------------| + 0 2 4 6 8 10 12 14 16 + */ + static piece force(int n, x_start x_start, x_stop x_stop, bin_trans bt) { + // TODO: Check that x_start < x_stop + const bool is_ordered = x_start.value() < x_stop.value(); + if (!is_ordered) + BOOST_THROW_EXCEPTION(std::invalid_argument("x_start must be < x_stop")); + + const width_start b0_obj = piece::calc_b0(x_start, x_stop, n, bt); + return piece(n, x_start, b0_obj, bt); } template T forward(T x) const noexcept { + static_assert(std::is_floating_point::value, "T must be a floating point type"); // Runs in hot loop, please measure impact of changes - switch (option_) { - case option::uniform: return forward_uniform(x); - case option::multiply: return forward_multiply(x); - case option::add: return forward_add(x); + T f_x; + if (bt_.is_uniform()) { + f_x = forward_uniform(x); + } else if (bt_.is_multiply()) { + f_x = forward_multiply(x); + } else if (bt_.is_add()) { + f_x = forward_add(x); + } else { + assert(false); // TODO: make unnecessary } + return f_x + accumulated_left_shift_; } - index_type size() const noexcept { return size_; } + index_type size() const noexcept { + return size_ic_ + accumulated_left_shift_ + accumulated_right_shift_; + } -private: - // Initialize - // Have: (x0, b0, n_bins) - // Calc: xN - - // Initialize - // Have: (x0, xN, n_bins) - // Calc: b0 - - // Add right - // Have: (x0, b0, n_bins++) - // Update: xN - - // Add left - // Have: (xN, b0, n_bins++) - // Calc: b0 - // Update: x0 - - void initialize(bin_stop stop) { - xN_ = x0_; - size_ = 0; - double b = b0_; - - while (!see_if_need_more_bins()) { - b = bin_trans_.calc_next_bin_size(b); - xN_ += b; - ++size_; + double x0() const noexcept { return x0_; } + double b0() const noexcept { return b0_; } + bin_trans bt() const noexcept { return bt_; } + index_type size_ic() const noexcept { return size_ic_; } + + internal_value_type xN() const noexcept { return xN_; } + + index_type accumulated_left_shift() const noexcept { return accumulated_left_shift_; } + index_type accumulated_right_shift() const noexcept { return accumulated_right_shift_; } + + // Bin ending at x_i has width + double calc_bin_width(double i) const noexcept { + const auto x_i = reverse(i); + const auto x_i_prev = reverse(i - 1); + return x_i - x_i_prev; + } + + static width_start calc_b0(x_start x0, x_stop xN, int N, bin_trans bt) { + double b0; + if (bt.is_uniform()) { + b0 = calc_b0_uniform(x0, xN, N); + } else if (bt.is_multiply()) { + b0 = calc_b0_multiply(x0, xN, N, bt); + } else if (bt.is_add()) { + b0 = calc_b0_add(x0, xN, N, bt); + } else { + assert(false); // TODO: make unnecessary } + return width_start(b0); } - bool see_if_need_more_bins() { - if (stop_.is_n_bins()) { - return size_ >= stop_.n_bins.value; - } else if (stop_.is_x()) { - return xN_ >= stop_.x.value; + template + internal_value_type reverse(T x) const noexcept { + if (bt_.is_uniform()) { + return reverse_uniform(x); + } else if (bt_.is_multiply()) { + return reverse_multiply(x); + } else if (bt_.is_add()) { + return reverse_add(x); + } else { + assert(false); // TODO: make unnecessary } } +private: + static piece rollout_right(int n, x_start xL, width_start b0, bin_trans bt) { + assert(0 < n); + return piece(n, xL, b0, bt); + } + + static piece rollout_left(int n_neg, x_start xR, width_start b0_R, + bin_trans bt_right) { + assert(n_neg < 0); + const int n_pos = -n_neg; + + // Roll out right + const piece piece_right(n_pos, xR, b0_R, bt_right); + + // Get width of last bin of right piece + const double width_last_bin = piece_right.calc_bin_width(n_pos); + + // Get location of last bin of right piece + const double x_last_bin = piece_right.reverse(n_pos); + + // Reflect last bin of right piece to get xL + const double delta = x_last_bin - xR.value(); + const double xL = xR.value() - delta; + + return piece(n_pos, x_start(xL), width_start(width_last_bin), bt_right.mirror()); + } + // Uniform spacing // // b0 b0 b0 b0 @@ -312,17 +308,30 @@ class piece { // x2 = x0 + b0 * 2 // x3 = x0 + b0 * 3 // ... - // x = x0 + b0 * N + // xN = x0 + b0 * N // - // Formula for N in terms of x - // N = (x - x0) / b0 + // Formula for N in terms of xN + // N = (xN - x0) / b0 // template - interval_value_type forward_linear(T x) const noexcept { + internal_value_type forward_uniform(T x) const noexcept { // Runs in hot loop, please measure impact of changes return ((x / unit_type{}) - x0_) / b0_; } + // Starting with: + // xN = x0 + b0 * N + // Solve for b0 to get: + // b0 = (xN - x0) / N + static double calc_b0_uniform(x_start x0, x_stop xN, int N) { + return (xN.value() - x0.value()) / N; + } + + template + internal_value_type reverse_uniform(T x) const noexcept { + return x0_ + x * b0_; + } + // Multiply spacing // // b0 b0*r b0*r² b0*r³ @@ -336,7 +345,7 @@ class piece { // x2 = r² * b0 / (r - 1) + ψ // x3 = r³ * b0 / (r - 1) + ψ // ... - // x = rᴺ * b0 / (r - 1) + ψ + // xN = rᴺ * b0 / (r - 1) + ψ // // Note: the first bin spacing is b0 // x1 - x0 = r * b0 / (r - 1) - 1 * b0 / (r - 1) @@ -344,20 +353,35 @@ class piece { // = b0 // // Find a formula for N - // x - x0 = b0 / (r - 1) * rᴺ - b0 / (r - 1) * 1 + // xN - x0 = b0 / (r - 1) * rᴺ - b0 / (r - 1) * 1 // = b0 / (r - 1) * (rᴺ - 1) // - // N in terms of x - // N = log2(((x - x0) / b0) * (r - 1) + 1) / log2(r) + // N in terms of xN + // N = log2(((xN - x0) / b0) * (r - 1) + 1) / log2(r) // template - interval_value_type forward_multiply(T x) const noexcept { + internal_value_type forward_multiply(T x) const noexcept { // Runs in hot loop, please measure impact of changes - const auto z = forward_linear(x); // z = (x - x0) / b0 - const auto r = bin_trans_.value(); + const auto z = forward_uniform(x); // z = (x - x0) / b0 + const auto r = bt_.value(); return std::log2(z * (r - 1) + 1) / std::log2(r); } + // Starting with: + // xN - x0 = b0 / (r - 1) * (rᴺ - 1) + // Solve for b0 to get: + // b0 = (xN - x0) * (r - 1) / (rᴺ - 1) + static double calc_b0_multiply(x_start x0, x_stop xN, int N, bin_trans bt) { + const auto r = bt.value(); + return (xN.value() - x0.value()) * (r - 1) / (std::pow(r, N) - 1); + } + + template + internal_value_type reverse_multiply(T x) const noexcept { + const auto r = bt_.value(); + return x0_ + b0_ * (std::pow(r, x) - 1) / (r - 1); + } + // Add spacing // // b0 b0+r b0+2r b0+3r @@ -375,35 +399,102 @@ class piece { // x = x0 + N * b0 + r * (1 + 2 + ... + N-1) // = x0 + N * b0 + r * (N * (N - 1) / 2) // - // N in terms of x - // 2 * (x - x0) / b0 = N * (N - 1) - // N² - N - 2 * (x - x0) / b0 = 0 + // Multiply by 2 + // 2 * x = 2 * x0 + 2 * N * b0 + r * (N * (N - 1)) + // Expand terms: + // 2 * x = 2 * x0 + 2 * N * b0 + rN² - rN + // Put all terms on one side: + // 0 = 2 * x0 + 2 * N * b0 + rN² - rN - 2 * x // Solve quadratic equation and take positive root + // N = (-1 + sqrt(1 + 8 * (x - x0) / b0 * r)) / 2 // template - interval_value_type forward_multiply(T x) const noexcept { + internal_value_type forward_add(T x) const noexcept { // Runs in hot loop, please measure impact of changes - const auto z = forward_linear(x); // z = (x - x0) / b0 - const auto r = bin_trans_.value(); - return (1 + std::sqrt(1 + 8 * z / r)) / 2; // Check copilot suggestion + const auto z = forward_uniform(x); // z = (x - x0) / b0 + const auto r = bt_.value(); + return (-1 + std::sqrt(1 + 8 * z * r)) / 2; + } + + // Starting with: + // 2 * x = 2 * x0 + 2 * N * b0 + r * (N * (N - 1)) + // Solve for b0 to get: + // b0 = 2 * (xN - x0) / (N * (N - 1)) - r + // + static double calc_b0_add(x_start x0, x_stop xN, int N, bin_trans bt) { + const auto r = bt.value(); + return 2 * (xN.value() - x0.value()) / (N * (N - 1)) - r; + } + + template + internal_value_type reverse_add(T x) const noexcept { + const auto r = bt_.value(); + return x0_ + x * b0_ + r * (x * (x - 1)) / 2; + } + + // Private constructor + piece(int n, x_start x0, width_start b0, bin_trans bt) + : size_ic_(n), x0_(x0.value()), b0_(b0.value()), bt_(bt) { + xN_ = reverse(size_ic_); } - internal_value_type x0_ { 0 } - double b0_{1}; - bin_trans bin_trans_{bin_trans::uniform{1}}; + index_type size_ic_{-9999}; + internal_value_type x0_{-9999.0}; + double b0_{-9999.0}; + bin_trans bt_{}; - internal_value_type xN_{0} index_type size_{0}; + internal_value_type xN_{-9999.0}; + + index_type accumulated_left_shift_{0}; + index_type accumulated_right_shift_{0}; }; -/// +/***Extrapolation and Attachment*** + + Pieces can be added to the right or left side of an existing piece. Pieces are + added with one of two approaches: extrapolation or attachment. Examples of each are + shown below. + Extrapolate 2 bins to the right side, doubling the bin size each time. + 0 2 4 6 8 + |---|---|---|---| + <- - - New piece - - -> + |---|---|---|---|-------|---------------| + 0 2 4 6 8 10 12 14 16 18 20 + Attach a new piece with a uniform bin spacing to the right side. + 0 2 4 6 8 + |---|---|---|---| + <- - - New piece - - -> + |---|---|---|---|-|-|-|-|-|-|-|-|-|-|-|-| + 0 2 4 6 8 10 12 14 16 18 20 + +***Syntax*** + +using P = piecewise; + +auto pa = P( // Create a new piece + P::x_start(1), // Start at x_start = 1 + P::width_start(1), // Start with a bin size of width_start = 1 + P::bin_trans::uniform(1), // Add bins of constant size 1 + P::stop::n_bins(4)); // Add 4 bins + +pa.add_right( // Add to the right side + P::bin_trans::multiply(2), // Extrapolate, doubling the bin size each time + P::stop::x(16)); // Stop when the piece contains x = 16 + +pa.add_left( // Add to the left side + P::width_start(1), // Start with a bin size of width_start = 1 + P::bin_trans::add(1), // Add bins of constant size 1 + P::stop::x(4)); // Stop when the piece contains x = 4 +*/ +template class piecewise { + using value_type = Value; + using unit_type = detail::get_unit_type; + using internal_value_type = detail::get_scale_type; + public: /// Construct a piecewise axis with an initial piece. - piecewise(x0 x0, b0 b0, bin_trans bin_trans, bin_stop stop) { - // - const auto p = piece(x0, b0, bin_trans, stop, true); - v_pieces_.push_back(p); - } + explicit piecewise(piece p) { v_pieces_.push_back(p); } /// template @@ -417,22 +508,6 @@ class piecewise { } } - /// - template - index_type index(ValueType x) const noexcept { - // Runs in hot loop, please measure impact of changes - const auto i = forward(x); - if (i < size_) { - if (0 <= i) - return static_cast(i); // 0 <= i < size - else - return -1; // i < 0 - } - // upper edge of last bin is inclusive if overflow bin is not present - if (!options_type::test(option::overflow) && z == 1) return size() - 1; - return size(); // also returned if x is NaN - } - /// Shifts the axis void shift_axis(index_type n) { if (n < 0) { @@ -443,78 +518,150 @@ class piecewise { } } - /// Returns index and shift (if axis has grown) for the passed argument. + void add_right(int n, bin_trans bt) { + double width_start = get_right_width_start(); + add_right(n, width_start, bt); + } + void add_right(int n, width_start width_start, bin_trans bt) { + double x_start = get_right_x(); + + const auto p = piece::rollout(n, x_start, width_start, bt); + v_pieces_.push_back(p); + } + + index_type size() const noexcept { + index_type size = 0; + for (const auto& p : v_pieces_) { size += p.size(); } + return size; + // One liner for the above function: + // return std::accumulate(v_pieces_.begin(), v_pieces_.end(), 0); + } + +private: + const piece& get_left() const noexcept { + assert(!v_pieces_.empty()); + return v_pieces_.front(); + } + + const piece& get_right() const noexcept { + assert(!v_pieces_.empty()); + return v_pieces_.back(); + } + + double get_left_x() const noexcept { return get_left().x0; } + double get_right_x() const noexcept { return get_right().x1; } + + double get_left_width_start() const noexcept { return get_left().width_start; } + double get_right_width_start() const noexcept { return get_right().width_start; } + + std::vector> v_pieces_{}; +}; + +/// +template +class int_shifter { + using value_type = Value; + +public: + /// Shifts the axis + void shift_axis(index_type n) { + if (n < 0) { + // Need to offset by this in the future + accumulated_left_shift_ += std::abs(n); + } else if (0 < n) { + // Need to offset by this in the future + accumulated_right_shift_ += n; + } + } + +private: + Payload payload_; + index_type accumulated_left_shift_{0}; + index_type accumulated_right_shift_{0}; +}; + +/// +template +class int_resolver { + using value_type = Value; + +public: + index_type index(value_type x) const noexcept { + // Runs in hot loop, please measure impact of changes + + y = payload_.index(x); + + if (y < size()) { + if (0 <= y) + return static_cast(y); // 0 <= i < size + else + return -1; // i < 0 + } + + // upper edge of last bin is inclusive if overflow bin is not present + if constexpr (std::is_floating_point::value) { + if (!options_type::test(option::overflow) && y == size()) return size() - 1; + } + + return size(); // also returned if x is NaN + } + std::pair update(value_type x) noexcept { - assert(options_type::test(option::growth)); - const auto i = forward(x); - if (i < size()) { - if (0 <= i) { - const auto i_int = static_cast(i); + // Runs in hot loop, please measure impact of changes + + y = payload_.index(x); + + if (y < size()) { + if (0 <= y) { + const auto i_int = static_cast(y); return {i_int, 0}; - } else - (i != -std::numeric_limits::infinity()) { - const auto i_int = static_cast(std::floor(i)); - shift_axis(i_int); - return {0, -i_int}; - } - else { // i is -infinity + } else if (y != -std::numeric_limits::infinity()) { + const auto i_int = static_cast(std::floor(y)); + shift_axis(i_int); + return {0, -i_int}; + } else { + // i is -infinity return {-1, 0}; } } // i either beyond range, infinite, or NaN - if (z < std::numeric_limits::infinity()) { - const auto i_int = static_cast(i); + if (y < std::numeric_limits::infinity()) { + const auto i_int = static_cast(y); const auto n = i_int - size() + 1; shift_axis(n); - return {i, -n}; + return {y, -n}; } // z either infinite or NaN return {size(), 0}; } - add_left(bin_trans bin_trans, bin_stop stop) { - double b0 = get_left_b0(); - add_left(b0, bin_trans, stop); - } - add_right(bin_trans bin_trans, bin_stop stop) { - double b0 = get_right_b0(); - add_right(b0, bin_trans, stop); - } - add_left(b0 b0, bin_trans bin_trans, bin_stop stop) { - double x0 = get_left_x(); - const auto p = piece(x0, b0, bin_trans, stop, false); - v_pieces_.insert(v_pieces_.begin(), p); - } - add_right(b0 b0, bin_trans bin_trans, bin_stop stop) { - double x0 = get_right_x(); - const auto p = piece(x0, b0, bin_trans, stop, true); - v_pieces_.push_back(p); - } + index_type size() const noexcept { return payload_.size(); } private: - double get_left_x() const { - assert(!v_pieces_.empty()); - return v_pieces_.front().x0; - } + Payload payload_; +}; - double get_right_x() const { - assert(!v_pieces_.empty()); - return v_pieces_.back().x1; - } +/// +template +class int_resolver_circular { + using value_type = Value; - /// Forward transform of external value x. - template - auto forward(T x) const { - return std::pow(x, power); - } +public: + index_type index(value_type x) const noexcept { + // Runs in hot loop, please measure impact of changes - /// Inverse transform of internal value x. - template - auto inverse(T x) const { - return std::pow(x, 1.0 / power); + y = payload_.index(x); + + // Not finite --> overflow bin + if constexpr (std::is_floating_point::value) { + if (!std::isfinite(y)) return payload_.size(); + } + + return y % payload_.size(); } - std::vector v_pieces_{}; +private: + Payload payload_; }; /** Axis for unit intervals on the real line. @@ -535,10 +682,9 @@ class piecewise { @tparam Options see boost::histogram::axis::option. */ template -class unit_regular - : public iterator_mixin>, - protected detail::replace_default, - public metadata_base_t { +class unit_regular : public iterator_mixin>, + // protected detail::replace_default, + public metadata_base_t { // these must be private, so that they are not automatically inherited using value_type = Value; using metadata_base = metadata_base_t; @@ -548,14 +694,16 @@ class unit_regular using unit_type = detail::get_unit_type; using internal_value_type = detail::get_scale_type; + using piecewise_type = piecewise; public: constexpr unit_regular() = default; - unit_regular(piecewise pw, metadata_type meta = {}, options_type options = {}) - : metadata_base(std::move(meta)), pw_(std::move(pw)) { + unit_regular(piecewise pw, metadata_type meta = {}, + options_type options = {}) + : metadata_base(std::move(meta)), piecewise_(std::move(pw)) { // static_asserts were moved here from class scope to satisfy deduction in gcc>=11 - static_assert(std::is_nothrow_move_constructible::value, + static_assert(std::is_nothrow_move_constructible::value, "piecewise must be no-throw move constructible"); static_assert(std::is_nothrow_move_assignable::value, "piecewise must be no-throw move assignable"); @@ -564,11 +712,6 @@ class unit_regular static_assert(!(options.test(option::circular) && options.test(option::growth)), "circular and growth options are mutually exclusive"); if (size() <= 0) BOOST_THROW_EXCEPTION(std::invalid_argument("bins > 0 required")); - if (!std::isfinite(min_) || !std::isfinite(delta_)) - BOOST_THROW_EXCEPTION( - std::invalid_argument("forward transform of start or stop invalid")); - if (delta_ == 0) - BOOST_THROW_EXCEPTION(std::invalid_argument("range of axis is zero")); } /// Constructor used by algorithm::reduce to shrink and rebin (not for users). @@ -629,7 +772,8 @@ class unit_regular private: index_type size_{0}; - internal_value_type min_{0} piecewise piecewise_{}; + internal_value_type min_{0}; + piecewise piecewise_{}; }; } // namespace axis diff --git a/include/boost/histogram/fwd.hpp b/include/boost/histogram/fwd.hpp index 3c48f784..3340840f 100644 --- a/include/boost/histogram/fwd.hpp +++ b/include/boost/histogram/fwd.hpp @@ -60,6 +60,12 @@ template class regular; +template +class piece; + +template +class unit_regular; + template class integer; diff --git a/test/CMakeLists.txt b/test/CMakeLists.txt index 337da698..cc4b0734 100644 --- a/test/CMakeLists.txt +++ b/test/CMakeLists.txt @@ -54,6 +54,7 @@ boost_test(TYPE run SOURCES axis_category_test.cpp) boost_test(TYPE run SOURCES axis_integer_test.cpp) boost_test(TYPE run SOURCES axis_option_test.cpp) boost_test(TYPE run SOURCES axis_regular_test.cpp) +boost_test(TYPE run SOURCES axis_piecewise_test.cpp) boost_test(TYPE run SOURCES axis_traits_test.cpp) boost_test(TYPE run SOURCES axis_variable_test.cpp) boost_test(TYPE run SOURCES axis_variant_test.cpp) diff --git a/test/Jamfile b/test/Jamfile index d38b5398..cef1a81c 100644 --- a/test/Jamfile +++ b/test/Jamfile @@ -57,6 +57,7 @@ alias cxx14 : [ run axis_integer_test.cpp ] [ run axis_option_test.cpp ] [ run axis_regular_test.cpp ] + [ run axis_piecewise_test.cpp ] [ run axis_traits_test.cpp ] [ run axis_variable_test.cpp ] [ run axis_variant_test.cpp ] diff --git a/test/axis_piecewise_test.cpp b/test/axis_piecewise_test.cpp new file mode 100644 index 00000000..c35f4455 --- /dev/null +++ b/test/axis_piecewise_test.cpp @@ -0,0 +1,392 @@ +// Copyright 2015-2019 Hans Dembinski +// +// Distributed under the Boost Software License, Version 1.0. +// (See accompanying file LICENSE_1_0.txt +// or copy at http://www.boost.org/LICENSE_1_0.txt) + +#include +#include +#include +#include +#include +#include +#include "axis.hpp" +#include "is_close.hpp" +#include "std_ostream.hpp" +#include "str.hpp" +#include "throw_exception.hpp" + +#include +#include +#include + +int main() { + using namespace boost::histogram; + // using def = use_default; + // namespace tr = axis::transform; + // namespace op = axis::option; + + // BOOST_TEST(std::is_nothrow_move_assignable>::value); + // BOOST_TEST(std::is_nothrow_move_constructible>::value); + + BOOST_TEST(std::is_nothrow_move_assignable>::value); + BOOST_TEST(std::is_nothrow_move_constructible>::value); + + // Construct rollout right + { + const auto p = axis::piece<>::rollout(5, axis::x_start(2.0), axis::width_start(1.1), + axis::bin_trans::multiply(1.2)); + BOOST_TEST_EQ(p.size(), 5); + BOOST_TEST(p.bt().is_multiply()); + BOOST_TEST_EQ(p.bt().value(), 1.2); + + BOOST_TEST_EQ(p.b0(), 1.1); + BOOST_TEST_EQ(p.x0(), 2.0); + BOOST_TEST_EQ(p.xN(), p.x0() + p.b0() * (1 + 1.2 + std::pow(1.2, 2) + + std::pow(1.2, 3) + std::pow(1.2, 4))); + + BOOST_TEST_EQ(p.size_ic(), 5); + BOOST_TEST_EQ(p.size(), 5); + BOOST_TEST_EQ(p.accumulated_left_shift(), 0); + BOOST_TEST_EQ(p.accumulated_right_shift(), 0); + + const auto b0_check = axis::piece<>::calc_b0(p.x0(), p.xN(), p.size_ic(), p.bt()); + BOOST_TEST_IS_CLOSE(b0_check.value(), p.b0(), 1e-12); + + BOOST_TEST_EQ(p.forward(2.0), 0.0); + BOOST_TEST_EQ(p.forward(2.0 + 1.1), 1.0); + BOOST_TEST_EQ(p.forward(2.0 + 1.1 * (1 + 1.2)), 2.0); + BOOST_TEST_EQ(p.forward(2.0 + 1.1 * (1 + 1.2 + std::pow(1.2, 2))), 3.0); + + for (int i = -100; i < 100; ++i) { + const auto d = i * 0.021; + BOOST_TEST_IS_CLOSE(p.reverse(p.forward(d)), d, 1e-12); + } + } + + // Construct rollout left + { + const auto p = axis::piece<>::rollout(-5, axis::x_start(2.0), axis::width_start(1.1), + axis::bin_trans::multiply(1.2)); + BOOST_TEST_EQ(p.size(), 5); + BOOST_TEST(p.bt().is_multiply()); + BOOST_TEST_EQ(p.bt().value(), 1 / 1.2); + + // const auto bt_rev = axis::bin_trans::multiply(1 / 1.2); + BOOST_TEST_IS_CLOSE(p.b0(), 1.1 * std::pow(1.2, 4), 1e-12); + BOOST_TEST_IS_CLOSE(p.xN(), 2.0, 1e-12); + + BOOST_TEST_IS_CLOSE( + p.x0(), + p.xN() - 1.1 * (1 + 1.2 + std::pow(1.2, 2) + std::pow(1.2, 3) + std::pow(1.2, 4)), + 1e-12); + + // const auto bt = p.bt(); + // BOOST_TEST(bt.is_multiply()); + // BOOST_TEST_EQ(bt.value(), 1.2); + + // BOOST_TEST_EQ(p.size_ic(), 5); + // BOOST_TEST_EQ(p.size(), 5); + // BOOST_TEST_EQ(p.accumulated_left_shift(), 0); + // BOOST_TEST_EQ(p.accumulated_right_shift(), 0); + + // std::cout << "p.bt().value(): " << p.bt().value() << std::endl; + // std::cout << "b0: " << p.b0() << std::endl; + // std::cout << "x0: " << p.x0() << std::endl; + // std::cout << "xN: " << p.xN() << std::endl; + + // for (int i = 1; i < 6; ++i) { + // std::cout << "bin_width: " << p.calc_bin_width(i) << std::endl; + // } + + // assert(false); + } + + // Construct piecewise + { + const auto p = axis::piece<>::rollout(-5, axis::x_start(2.0), axis::width_start(1.1), + axis::bin_trans::multiply(1.2)); + auto pw = axis::piecewise(p); + } + + // // bad_ctors + // { + // BOOST_TEST_THROWS(axis::unit_regular<>(1, 0, 0), std::invalid_argument); + // BOOST_TEST_THROWS(axis::unit_regular<>(0, 0, 1), std::invalid_argument); + // } + + // // ctors and assignment + // { + // axis::unit_regular<> a{4, -2, 2}; + // axis::unit_regular<> b; + // BOOST_TEST_NE(a, b); + // b = a; + // BOOST_TEST_EQ(a, b); + // axis::unit_regular<> c = std::move(b); + // BOOST_TEST_EQ(c, a); + // axis::unit_regular<> d; + // BOOST_TEST_NE(c, d); + // d = std::move(c); + // BOOST_TEST_EQ(d, a); + // } + + // // input, output + // { + // axis::unit_regular<> a{4, -2, 2, "foo"}; + // BOOST_TEST_EQ(a.metadata(), "foo"); + // const auto& cref = a; + // BOOST_TEST_EQ(cref.metadata(), "foo"); + // cref.metadata() = "bar"; // this is allowed + // BOOST_TEST_EQ(cref.metadata(), "bar"); + // BOOST_TEST_EQ(a.value(0), -2); + // BOOST_TEST_EQ(a.value(1), -1); + // BOOST_TEST_EQ(a.value(2), 0); + // BOOST_TEST_EQ(a.value(3), 1); + // BOOST_TEST_EQ(a.value(4), 2); + // BOOST_TEST_EQ(a.bin(-1).lower(), -std::numeric_limits::infinity()); + // BOOST_TEST_EQ(a.bin(-1).upper(), -2); + // BOOST_TEST_EQ(a.bin(a.size()).lower(), 2); + // BOOST_TEST_EQ(a.bin(a.size()).upper(), std::numeric_limits::infinity()); + // BOOST_TEST_EQ(a.index(-10.), -1); + // BOOST_TEST_EQ(a.index(-2.1), -1); + // BOOST_TEST_EQ(a.index(-2.0), 0); + // BOOST_TEST_EQ(a.index(-1.1), 0); + // BOOST_TEST_EQ(a.index(0.0), 2); + // BOOST_TEST_EQ(a.index(0.9), 2); + // BOOST_TEST_EQ(a.index(1.0), 3); + // BOOST_TEST_EQ(a.index(10.), 4); + // BOOST_TEST_EQ(a.index(-std::numeric_limits::infinity()), -1); + // BOOST_TEST_EQ(a.index(std::numeric_limits::infinity()), 4); + // BOOST_TEST_EQ(a.index(std::numeric_limits::quiet_NaN()), 4); + + // BOOST_TEST_EQ(str(a), + // "regular(4, -2, 2, metadata=\"bar\", options=underflow | overflow)"); + // } + + // // with inverted range + // { + // axis::unit_regular<> a{2, 1, -2}; + // BOOST_TEST_EQ(a.bin(-1).lower(), std::numeric_limits::infinity()); + // BOOST_TEST_EQ(a.bin(0).lower(), 1); + // BOOST_TEST_EQ(a.bin(1).lower(), -0.5); + // BOOST_TEST_EQ(a.bin(2).lower(), -2); + // BOOST_TEST_EQ(a.bin(2).upper(), -std::numeric_limits::infinity()); + // BOOST_TEST_EQ(a.index(2), -1); + // BOOST_TEST_EQ(a.index(1.001), -1); + // BOOST_TEST_EQ(a.index(1), 0); + // BOOST_TEST_EQ(a.index(0), 0); + // BOOST_TEST_EQ(a.index(-0.499), 0); + // BOOST_TEST_EQ(a.index(-0.5), 1); + // BOOST_TEST_EQ(a.index(-1), 1); + // BOOST_TEST_EQ(a.index(-2), 2); + // BOOST_TEST_EQ(a.index(-20), 2); + // } + + // // with log transform + // { + // auto a = axis::unit_regular{2, 1e0, 1e2}; + // BOOST_TEST_EQ(a.bin(-1).lower(), 0.0); + // BOOST_TEST_IS_CLOSE(a.bin(0).lower(), 1.0, 1e-9); + // BOOST_TEST_IS_CLOSE(a.bin(1).lower(), 10.0, 1e-9); + // BOOST_TEST_IS_CLOSE(a.bin(2).lower(), 100.0, 1e-9); + // BOOST_TEST_EQ(a.bin(2).upper(), std::numeric_limits::infinity()); + + // BOOST_TEST_EQ(a.index(-1), 2); // produces NaN in conversion + // BOOST_TEST_EQ(a.index(0), -1); + // BOOST_TEST_EQ(a.index(1), 0); + // BOOST_TEST_EQ(a.index(9), 0); + // BOOST_TEST_EQ(a.index(10), 1); + // BOOST_TEST_EQ(a.index(90), 1); + // BOOST_TEST_EQ(a.index(100), 2); + // BOOST_TEST_EQ(a.index(std::numeric_limits::infinity()), 2); + + // BOOST_TEST_THROWS((axis::unit_regular{2, -1, 0}), + // std::invalid_argument); + + // BOOST_TEST_CSTR_EQ( + // str(a).c_str(), + // "regular(transform::log{}, 2, 1, 100, options=underflow | overflow)"); + // } + + // // with sqrt transform + // { + // axis::unit_regular a(2, 0, 4); + // // this is weird, but -inf * -inf = inf, thus the lower bound + // BOOST_TEST_EQ(a.bin(-1).lower(), std::numeric_limits::infinity()); + // BOOST_TEST_IS_CLOSE(a.bin(0).lower(), 0.0, 1e-9); + // BOOST_TEST_IS_CLOSE(a.bin(1).lower(), 1.0, 1e-9); + // BOOST_TEST_IS_CLOSE(a.bin(2).lower(), 4.0, 1e-9); + // BOOST_TEST_EQ(a.bin(2).upper(), std::numeric_limits::infinity()); + + // BOOST_TEST_EQ(a.index(-1), 2); // produces NaN in conversion + // BOOST_TEST_EQ(a.index(0), 0); + // BOOST_TEST_EQ(a.index(0.99), 0); + // BOOST_TEST_EQ(a.index(1), 1); + // BOOST_TEST_EQ(a.index(3.99), 1); + // BOOST_TEST_EQ(a.index(4), 2); + // BOOST_TEST_EQ(a.index(100), 2); + // BOOST_TEST_EQ(a.index(std::numeric_limits::infinity()), 2); + + // BOOST_TEST_EQ(str(a), + // "regular(transform::sqrt{}, 2, 0, 4, options=underflow | overflow)"); + // } + + // // with pow transform + // { + // axis::unit_regular a(tr::pow{0.5}, 2, 0, 4); + // // this is weird, but -inf * -inf = inf, thus the lower bound + // BOOST_TEST_EQ(a.bin(-1).lower(), std::numeric_limits::infinity()); + // BOOST_TEST_IS_CLOSE(a.bin(0).lower(), 0.0, 1e-9); + // BOOST_TEST_IS_CLOSE(a.bin(1).lower(), 1.0, 1e-9); + // BOOST_TEST_IS_CLOSE(a.bin(2).lower(), 4.0, 1e-9); + // BOOST_TEST_EQ(a.bin(2).upper(), std::numeric_limits::infinity()); + + // BOOST_TEST_EQ(a.index(-1), 2); // produces NaN in conversion + // BOOST_TEST_EQ(a.index(0), 0); + // BOOST_TEST_EQ(a.index(0.99), 0); + // BOOST_TEST_EQ(a.index(1), 1); + // BOOST_TEST_EQ(a.index(3.99), 1); + // BOOST_TEST_EQ(a.index(4), 2); + // BOOST_TEST_EQ(a.index(100), 2); + // BOOST_TEST_EQ(a.index(std::numeric_limits::infinity()), 2); + + // BOOST_TEST_EQ(str(a), + // "regular(transform::pow{0.5}, 2, 0, 4, options=underflow | + // overflow)"); + // } + + // // with step + // { + // axis::unit_regular<> a(axis::step(0.5), 1, 3); + // BOOST_TEST_EQ(a.size(), 4); + // BOOST_TEST_EQ(a.bin(-1).lower(), -std::numeric_limits::infinity()); + // BOOST_TEST_EQ(a.value(0), 1); + // BOOST_TEST_EQ(a.value(1), 1.5); + // BOOST_TEST_EQ(a.value(2), 2); + // BOOST_TEST_EQ(a.value(3), 2.5); + // BOOST_TEST_EQ(a.value(4), 3); + // BOOST_TEST_EQ(a.bin(4).upper(), std::numeric_limits::infinity()); + + // axis::unit_regular<> b(axis::step(0.5), 1, 3.1); + // BOOST_TEST_EQ(a, b); + // } + + // // with circular option + // { + // axis::circular<> a{4, 0, 1}; + // BOOST_TEST_EQ(a.bin(-1).lower(), a.bin(a.size() - 1).lower() - 1); + // BOOST_TEST_EQ(a.index(-1.0 * 3), 0); + // BOOST_TEST_EQ(a.index(0.0), 0); + // BOOST_TEST_EQ(a.index(0.25), 1); + // BOOST_TEST_EQ(a.index(0.5), 2); + // BOOST_TEST_EQ(a.index(0.75), 3); + // BOOST_TEST_EQ(a.index(1.0), 0); + // BOOST_TEST_EQ(a.index(std::numeric_limits::infinity()), 4); + // BOOST_TEST_EQ(a.index(-std::numeric_limits::infinity()), 4); + // BOOST_TEST_EQ(a.index(std::numeric_limits::quiet_NaN()), 4); + // } + + // // with growth + // { + // using pii_t = std::pair; + // axis::unit_regular a{1, 0, 1}; + // BOOST_TEST_EQ(a.size(), 1); + // BOOST_TEST_EQ(a.update(0), pii_t(0, 0)); + // BOOST_TEST_EQ(a.size(), 1); + // BOOST_TEST_EQ(a.update(1), pii_t(1, -1)); + // BOOST_TEST_EQ(a.size(), 2); + // BOOST_TEST_EQ(a.value(0), 0); + // BOOST_TEST_EQ(a.value(2), 2); + // BOOST_TEST_EQ(a.update(-1), pii_t(0, 1)); + // BOOST_TEST_EQ(a.size(), 3); + // BOOST_TEST_EQ(a.value(0), -1); + // BOOST_TEST_EQ(a.value(3), 2); + // BOOST_TEST_EQ(a.update(-10), pii_t(0, 9)); + // BOOST_TEST_EQ(a.size(), 12); + // BOOST_TEST_EQ(a.value(0), -10); + // BOOST_TEST_EQ(a.value(12), 2); + // BOOST_TEST_EQ(a.update(std::numeric_limits::infinity()), pii_t(a.size(), + // 0)); BOOST_TEST_EQ(a.update(std::numeric_limits::quiet_NaN()), + // pii_t(a.size(), 0)); + // BOOST_TEST_EQ(a.update(-std::numeric_limits::infinity()), pii_t(-1, 0)); + // } + + // // axis with overflow bin represents open interval + // { + // axis::unit_regular a{2, 0, 1}; + // BOOST_TEST_EQ(a.index(0), 0); + // BOOST_TEST_EQ(a.index(0.49), 0); + // BOOST_TEST_EQ(a.index(0.50), 1); + // BOOST_TEST_EQ(a.index(0.99), 1); + // BOOST_TEST_EQ(a.index(1), 2); // overflow bin + // BOOST_TEST_EQ(a.index(1.1), 2); // overflow bin + // } + + // // axis without overflow bin represents a closed interval + // { + // axis::unit_regular a{2, 0, 1}; + // BOOST_TEST_EQ(a.index(0), 0); + // BOOST_TEST_EQ(a.index(0.49), 0); + // BOOST_TEST_EQ(a.index(0.50), 1); + // BOOST_TEST_EQ(a.index(0.99), 1); + // BOOST_TEST_EQ(a.index(1), 1); // last ordinary bin + // BOOST_TEST_EQ(a.index(1.1), 2); // out of range + // } + + // // iterators + // { + // test_axis_iterator(axis::unit_regular<>(5, 0, 1), 0, 5); + // test_axis_iterator(axis::unit_regular(5, 0, 1), 0, + // 5); test_axis_iterator(axis::circular<>(5, 0, 1), 0, 5); + // } + + // // bin_type streamable + // { + // auto test = [](const auto& x, const char* ref) { + // std::ostringstream os; + // os << x; + // BOOST_TEST_EQ(os.str(), std::string(ref)); + // }; + + // auto a = axis::unit_regular<>(2, 0, 1); + // test(a.bin(0), "[0, 0.5)"); + // } + + // // null_type streamable + // { + // auto a = axis::unit_regular(2, 0, 1); + // BOOST_TEST_EQ(str(a), "regular(2, 0, 1, options=underflow | overflow)"); + // } + + // // shrink and rebin + // { + // using A = axis::unit_regular<>; + // auto a = A(5, 0, 5); + // auto b = A(a, 1, 4, 1); + // BOOST_TEST_EQ(b.size(), 3); + // BOOST_TEST_EQ(b.value(0), 1); + // BOOST_TEST_EQ(b.value(3), 4); + // auto c = A(a, 0, 4, 2); + // BOOST_TEST_EQ(c.size(), 2); + // BOOST_TEST_EQ(c.value(0), 0); + // BOOST_TEST_EQ(c.value(2), 4); + // auto e = A(a, 1, 5, 2); + // BOOST_TEST_EQ(e.size(), 2); + // BOOST_TEST_EQ(e.value(0), 1); + // BOOST_TEST_EQ(e.value(2), 5); + // } + + // // shrink and rebin with circular option + // { + // using A = axis::circular<>; + // auto a = A(4, 1, 5); + // BOOST_TEST_THROWS(A(a, 1, 4, 1), std::invalid_argument); + // BOOST_TEST_THROWS(A(a, 0, 3, 1), std::invalid_argument); + // auto b = A(a, 0, 4, 2); + // BOOST_TEST_EQ(b.size(), 2); + // BOOST_TEST_EQ(b.value(0), 1); + // BOOST_TEST_EQ(b.value(2), 5); + // } + + return boost::report_errors(); +} From 821ea49d9d898ca49eb63be0c56a2dad865b1af7 Mon Sep 17 00:00:00 2001 From: Ryan Date: Tue, 13 Jun 2023 15:04:58 -0400 Subject: [PATCH 03/10] divided into piece, piecewise, + int_resolver --- include/boost/histogram/axis.hpp | 2 + include/boost/histogram/axis/int_resolver.hpp | 214 +++++ include/boost/histogram/axis/piece.hpp | 512 +++++++++++ include/boost/histogram/axis/piecewise.hpp | 818 +++--------------- test/CMakeLists.txt | 2 + test/Jamfile | 2 + test/axis_int_resolver_test.cpp | 111 +++ test/axis_piece_test.cpp | 213 +++++ test/axis_piecewise_test.cpp | 401 +-------- 9 files changed, 1243 insertions(+), 1032 deletions(-) create mode 100644 include/boost/histogram/axis/int_resolver.hpp create mode 100644 include/boost/histogram/axis/piece.hpp create mode 100644 test/axis_int_resolver_test.cpp create mode 100644 test/axis_piece_test.cpp diff --git a/include/boost/histogram/axis.hpp b/include/boost/histogram/axis.hpp index 078be5cf..bc95f80c 100644 --- a/include/boost/histogram/axis.hpp +++ b/include/boost/histogram/axis.hpp @@ -19,7 +19,9 @@ #include #include +#include #include +#include #include #include #include diff --git a/include/boost/histogram/axis/int_resolver.hpp b/include/boost/histogram/axis/int_resolver.hpp new file mode 100644 index 00000000..a9b7e130 --- /dev/null +++ b/include/boost/histogram/axis/int_resolver.hpp @@ -0,0 +1,214 @@ +// Copyright 2015-2018 Hans Dembinski +// +// Distributed under the Boost Software License, Version 1.0. +// (See accompanying file LICENSE_1_0.txt +// or copy at http://www.boost.org/LICENSE_1_0.txt) + +#ifndef BOOST_HISTOGRAM_AXIS_INT_RESOLVER_HPP +#define BOOST_HISTOGRAM_AXIS_INT_RESOLVER_HPP + +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include + +namespace boost { +namespace histogram { +namespace axis { + +/** Bin shift integrator + + Accumulates bin shifts to support growable axes. +*/ +template +class bin_shift_integrator { + using value_type = Value; + using internal_value_type = detail::get_scale_type; + +public: + /// Constructor for bin_shift_integrator + explicit bin_shift_integrator(const Payload& payload) : payload_(payload) {} + + /// Shifts the axis + void shift_axis(index_type n) { + if (n < 0) { + bins_under_ += std::abs(n); + } else if (0 < n) { + bins_over_ += n; + } + } + + /// The mapping from bin space Y to input space X + template + T inverse(T y) const noexcept { + return payload_.inverse(y - bins_under_); + } + + /// The mapping from input space X to bin space Y + template + T forward(T x) const noexcept { + return payload_.forward(x) + bins_under_; + } + + /// The number of bins in the axis + index_type size() const noexcept { return payload_.size() + bins_under_ + bins_over_; } + +private: + Payload payload_; + index_type bins_under_{0}; + index_type bins_over_{0}; +}; + +/** Int resolver + + Resolves float bin numbers to integer bin numbers. +*/ +template +class int_resolver_linear { + using value_type = Value; + using unit_type = detail::get_unit_type; + using internal_value_type = detail::get_scale_type; + +public: + /// Constructor for int_resolver_linear + explicit int_resolver_linear(const Payload& payload) : payload_(payload) {} + + /// The mapping from input space X to integer bin space Y + index_type index(value_type x) const noexcept { + // Runs in hot loop, please measure impact of changes + + const value_type y = payload_.forward(x); + + if (y < size()) { + if (0 <= y) + return static_cast(y); // 0 <= i < size + else + return -1; // i < 0 + } + + // upper edge of last bin is inclusive if overflow bin is not present + if (std::is_floating_point::value) { + // TODO: enable support for this feature + // if (!options_type::test(option::overflow) && y == size()) return size() - 1; + } + + return size(); // also returned if x is NaN + } + + /// TODO: move this to a new class. + std::pair update(value_type x) noexcept { + const value_type y = forward(x); + + if (y < size()) { + if (0 <= y) { + const auto i_int = static_cast(y); + return {i_int, 0}; + } else if (y != -std::numeric_limits::infinity()) { + const auto i_int = static_cast(std::floor(y)); + payload_.shift_axis(i_int); + return {0, -i_int}; + } else { + return {-1, 0}; // i is -infinity + } + } + // i either beyond range, infinite, or NaN + if (y < std::numeric_limits::infinity()) { + const auto i_int = static_cast(y); + const auto n = i_int - size() + 1; + payload_.shift_axis(n); + return {y, -n}; + } + // z either infinite or NaN + return {size(), 0}; + } + + /// The mapping from bin space Y to input space X + value_type value(real_index_type i) const noexcept { + if (i < 0) { + return -std::numeric_limits::infinity(); + } else if (i <= size()) { + return static_cast(payload_.inverse(i) * unit_type{}); + } else { + return std::numeric_limits::infinity(); + } + } + + index_type size() const noexcept { return payload_.size(); } + +private: + internal_value_type forward(internal_value_type x) const noexcept { + return payload_.forward(x); + } + + Payload payload_; +}; + +/** Circular int resolver + + Resolves float bin numbers to integer bin numbers. +*/ +template +class int_resolver_circular { + using value_type = Value; + +public: + /// Constructor for int_resolver_circular + explicit int_resolver_circular(const Payload& payload, Value x_low, Value x_high) + : payload_(payload) + , x_low_(x_low) + , x_high_(x_high) + , x_low_plus_delta_x_(x_high - x_low) { + if (x_high <= x_low) BOOST_THROW_EXCEPTION(std::invalid_argument("x_high <= x_low")); + if (payload.x0() < x_low) + BOOST_THROW_EXCEPTION(std::invalid_argument("payload.x0() < x_low")); + if (x_high < payload.xN()) + BOOST_THROW_EXCEPTION(std::invalid_argument("x_high < payload.xN()")); + } + + index_type index(value_type x) const noexcept { + // Take the mod of the input + x -= x_low_plus_delta_x_ * std::floor((x - x_low_) / x_low_plus_delta_x_); + + if (x < x_low_ || x_high_ < x) { + return payload_.size(); + } else { + const value_type y = payload_.forward(x); + + if (std::isfinite(y)) { return static_cast(y); } + + return payload_.size(); + } + } + +private: + Payload payload_; + Value x_low_; + Value x_high_; + Value x_low_plus_delta_x_; +}; + +} // namespace axis +} // namespace histogram +} // namespace boost + +#endif diff --git a/include/boost/histogram/axis/piece.hpp b/include/boost/histogram/axis/piece.hpp new file mode 100644 index 00000000..92098bb9 --- /dev/null +++ b/include/boost/histogram/axis/piece.hpp @@ -0,0 +1,512 @@ +// Copyright 2015-2018 Hans Dembinski +// +// Distributed under the Boost Software License, Version 1.0. +// (See accompanying file LICENSE_1_0.txt +// or copy at http://www.boost.org/LICENSE_1_0.txt) + +#ifndef BOOST_HISTOGRAM_AXIS_PIECE_HPP +#define BOOST_HISTOGRAM_AXIS_PIECE_HPP + +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include + +namespace boost { +namespace histogram { +namespace axis { + +/** The abstract piece class that other pieces inherit from + + The class has two pure virtual functions: + 1. inverse() maps from bin space Y to input space X + 2. forward() maps from input space X to bin space Y + + This abstract class has three data members: + 1. size_ic_ number of bins + 2. x0_ lower edge of the first bin + 3. xN_ upper edge of the last bin + + The member x0_ is used internally by the forward and inverse functions. The member + xN_ can be used by other classes to avoid calling inverse() to get the location of + on the last bin. +*/ +template +class piece_base { + using value_type = Value; + using internal_value_type = detail::get_scale_type; + +public: + /// The mapping from bin space Y to input space X + virtual internal_value_type inverse(internal_value_type y) const noexcept = 0; + + /// The mapping from input space X to bin space Y + virtual internal_value_type forward(internal_value_type x) const noexcept = 0; + + /// Creates a piece for the given arguments. This function checks the orientation of + /// the transform and makes sure that transforms are not inverted. + template + static PieceType create(Args... args) { + auto p = PieceType(args...); + p.xN_ = p.inverse(p.size_ic_); // Set xN_ + + assert(0 < p.size_ic_); + assert(p.x0_ < p.xN_); + return p; + } + + /// Calculates the width of bins. Calling calc_bin_width(0) gives the width of the first + /// bin. + double calc_bin_width(double i) const noexcept { + const auto x_i_next = inverse(i + 1); + const auto x_i = inverse(i); + return x_i_next - x_i; + } + + /// Calculates the width of the last bin. + double calc_bin_width_last() const noexcept { return calc_bin_width(size_ic_ - 1); } + + /// The number of bins + index_type size() const noexcept { return size_ic_; } + + /// The lower edge of the first bin (i.e., bin 0). + internal_value_type x0() const noexcept { return x0_; } + + /// The upper edge of the last bin (i.e., bin N - 1). + internal_value_type xN() const noexcept { return xN_; } + +protected: + // Constructor called by derived classes + piece_base(int N, double x0) : size_ic_{N}, x0_{x0} {} + +private: + index_type size_ic_{-9999}; + internal_value_type x0_{-9999.0}; + internal_value_type xN_{-9999.0}; +}; + +/** Piece for a unit transformation. + +*/ +template +class piece_unit : public piece_base> { + using value_type = Value; + using internal_value_type = detail::get_scale_type; + +public: + /// The mapping from bin space Y to input space X + internal_value_type inverse(internal_value_type y) const noexcept override { + return y + this->x0(); + } + + /// The mapping from input space X to bin space Y + internal_value_type forward(internal_value_type x) const noexcept override { + return x - this->x0(); + } + +private: + // Private constructor below is called by the base class. + friend class piece_base>; + + piece_unit(int N, double x0) : piece_base>(N, x0) {} +}; + +/** Piece for a variable transformation. + +*/ +template +class piece_variable : public piece_base> { + using value_type = Value; + using internal_value_type = detail::get_scale_type; + +public: + /// The mapping from bin space Y to input space X + internal_value_type inverse(internal_value_type y) const noexcept override { + // NOTE: copied from variable.hpp + if (y < 0) return detail::lowest(); + if (y == this->size()) return vec_.back(); + if (y > this->size()) return detail::highest(); + const auto k = static_cast(y); // precond: y >= 0 + const real_index_type z = y - k; + // check z == 0 needed to avoid returning nan when vec_[k + 1] is infinity + return (1.0 - z) * vec_[k] + (z == 0 ? 0 : z * vec_[k + 1]); + } + + /// The mapping from input space X to bin space Y + internal_value_type forward(internal_value_type x) const noexcept override { + // NOTE: copied from variable.hpp + if (x == vec_.back()) return this->size() - 1; // TODO: support growth option + return static_cast(std::upper_bound(vec_.begin(), vec_.end(), x) - + vec_.begin() - 1); + } + +private: + // Private constructor below is called by the base class. + friend class piece_base>; + + piece_variable(const std::vector& vec) + : piece_base>(vec.size() - 1, vec[0]), vec_{vec} {} + + std::vector vec_; +}; + +/** Abstract variable width piece class + +*/ +template +class piece_b0 : public piece_base { + using value_type = Value; + using internal_value_type = detail::get_scale_type; + +public: + /** Finds the piece whose size is known, but bin spacing is not. + + @param n number of bins + @param x0 starting point + @param xN stopping point + @param args additional arguments + + The constructor throws `std::invalid_argument` if: + 1. the bin number is zero, or + 2. x0 <= xN. + */ + template + static PieceType solve_b0(int n, double x0, double xN, Args... args) { + // TODO: Check that x0 < xN + const bool is_ordered = x0 < xN; + if (!is_ordered) BOOST_THROW_EXCEPTION(std::invalid_argument("x0 must be < xN")); + + const double b0_obj = PieceType::calc_b0(n, x0, xN, args...); + return PieceType::create(n, x0, b0_obj, args...); + } + + /// Returns the width of the first bin. + double b0() const noexcept { return b0_; } + + piece_b0(int N, double x0, double b0) : piece_base(N, x0), b0_(b0) {} + +private: + friend class piece_base; + + double b0_; +}; + +/** Piece for uniform spacing + + Uniform bin spacing b0 + + b0 b0 b0 b0 + y=0 y=1 y=2 y=3 y=4 + |------|------|------|------| + x0 x4 + + The sequence + x0 = x0 + x1 = x0 + b0 + x2 = x0 + b0 * 2 + x3 = x0 + b0 * 3 + ... + xN = x0 + b0 * yN + +*/ +template +class piece_uniform : public piece_b0> { + using value_type = Value; + using unit_type = detail::get_unit_type; + using internal_value_type = detail::get_scale_type; + +public: + /** The mapping from bin space Y to input space X. + + Uses the formula: + x = x0 + b0 * y + */ + // Solves for x in terms of y + internal_value_type inverse(internal_value_type y) const noexcept override { + return this->x0() + y * this->b0(); + } + + /** The mapping from input space X to bin space Y. + + Solving x = x0 + b0 * y for y gives the formula: + y = (x - x0) / b0 + */ + internal_value_type forward(internal_value_type x) const noexcept override { + return ((x / unit_type{}) - this->x0()) / this->b0(); + } + + /** Calculates the bin width. + + Uses the formula: + b0 = (xN - x0) / N + */ + static double calc_b0(int N, double x0, double xN) { return (xN - x0) / N; } + + /// Calculates the next bin width. The bin width is constant. + static double next_b(double b0) { return b0; } + +private: + // Private constructor below is called by the base class. + friend class piece_base>; + + piece_uniform(int N, double x0, double b0) + : piece_b0>(N, x0, b0) {} +}; + +/** Creates a piece where the bin size is multiplied by a constant + + b0 b0*r b0*r² b0*r³ + N=0 N=1 N=2 N=3 N=4 + |---|-----|-----------|------------------| + x0 x4 + + The sequence (for some value of ψ) + x0 = 1 * b0 / (r - 1) + ψ + x1 = r * b0 / (r - 1) + ψ + x2 = r² * b0 / (r - 1) + ψ + x3 = r³ * b0 / (r - 1) + ψ + ... + xN = rᴺ * b0 / (r - 1) + ψ + + Note: the first bin spacing is b0 + x1 - x0 = r * b0 / (r - 1) - 1 * b0 / (r - 1) + = (r - 1) * b0 / (r - 1) + = b0 + + Find a formula for N + xN - x0 = b0 / (r - 1) * rᴺ - b0 / (r - 1) * 1 + = b0 / (r - 1) * (rᴺ - 1) + + N in terms of xN + N = log2(((xN - x0) / b0) * (r - 1) + 1) / log2(r) + +*/ +template +class piece_multiply : public piece_b0> { + using value_type = Value; + using unit_type = detail::get_unit_type; + using internal_value_type = detail::get_scale_type; + +public: + /** The mapping from bin space Y to input space X. + + Uses the formula: + x = x0 + b0 * (rᴺ - 1) / (r - 1) + */ + internal_value_type inverse(internal_value_type N) const noexcept override { + return this->x0() + this->b0() * (std::pow(r_, N) - 1) / (r_ - 1); + } + + /** The mapping from input space X to bin space Y. + + Solving the inverse formula for N gives the formula: + N = log2(((x - x0) / b0) * (r - 1) + 1) / log2(r) + */ + internal_value_type forward(internal_value_type x) const noexcept override { + const auto z = ((x / unit_type{}) - this->x0()) / this->b0(); + return std::log2(z * (r_ - 1) + 1) / std::log2(r_); + } + + /** Calculates the bin width for the first bin. + + Uses the formula: + b0 = (xN - x0) * (r - 1) / (rᴺ - 1) + */ + static double calc_b0(int N, double x0, double xN, double r) { + return (xN - x0) * (r - 1) / (std::pow(r, N) - 1); + } + + /// Calculates the next bin width. + static double next_b(double b0, double r) { return b0 * r; } + + /// Returns the ratio between bin widths. + double r() const noexcept { return r_; } + +private: + // Private constructor below is called by the base class. + friend class piece_base>; + + piece_multiply(int N, double x0, double b0, double r) + : piece_b0>(N, x0, b0), r_(r) {} + + double r_; +}; + +/** Piece for adding a constant to the bin spacing + + b0 b0+r b0+2r b0+3r + N=0 N=1 N=2 N=3 N=4 + |---|-----|-----------|------------------| + x0 x4 + + The sequence + x0 = x0 = 0 * b0 + x1 = x0 + 1 * b0 + r * 0 + x2 = x0 + 2 * b0 + r * 1 + x3 = x0 + 3 * b0 + r * (1 + 2) + x4 = x0 + 4 * b0 + r * (1 + 2 + 3) + ... + x = x0 + N * b0 + r * (1 + 2 + ... + N-1) + x = x0 + N * b0 + r * (N * (N - 1) / 2) + +*/ +template +class piece_add : public piece_b0> { + using value_type = Value; + using internal_value_type = detail::get_scale_type; + +public: + /// The mapping from bin space Y to input space X + internal_value_type inverse(internal_value_type N) const noexcept override { + return this->x0() + N * this->b0() + r_ * (N * (N - 1)) / 2; + } + + /** The mapping from input space X to bin space Y. + + Starting with: + x = x0 + N * b0 + r * (N * (N - 1) / 2) + Multiply by 2: + 2 * x = 2 * x0 + 2 * N * b0 + r * (N * (N - 1)) + Expand terms: + 2 * x = 2 * x0 + 2 * N * b0 + rN² - rN + Collect terms by powers of N: + 0 = rN² + (2 * b0 - r)N + (2 * x0 - 2 * x) + + a = r + b = 2 * b0 - r + c = 2 * x0 - 2 * x + */ + internal_value_type forward(internal_value_type x) const noexcept override { + const double a = r_; + const double b = 2 * this->b0() - r_; + const double c = 2 * this->x0() - 2 * x; + const auto roots = quadratic_roots(a, b, c); + return roots.second; // Take the positive root + } + + /// Finds the roots of the quadratic equation ax² + bx + c = 0. + static std::pair quadratic_roots(const value_type& a, + const value_type& b, + const value_type& c) noexcept { + // https://people.csail.mit.edu/bkph/articles/Quadratics.pdf + + const value_type two_a = 2 * a; + const value_type two_c = 2 * c; + const value_type sqrt_bb_4ac = std::sqrt(b * b - two_a * two_c); + + if (b >= 0) { + const value_type root1 = (-b - sqrt_bb_4ac) / two_a; + const value_type root2 = two_c / (-b - sqrt_bb_4ac); + return {root1, root2}; + } else { + const value_type root1 = two_c / (-b + sqrt_bb_4ac); + const value_type root2 = (-b + sqrt_bb_4ac) / two_a; + return {root1, root2}; + } + } + + /** Calculates the bin width for the first bin. + + Uses the formula: + b0 = (xN - x0 - r * (N * (N - 1) / 2)) / N + */ + static double calc_b0(int N, double x0, double xN, double r) { + return (xN - x0 - r * (N * (N - 1) / 2)) / N; + } + + /// Calculates the next bin width. + static double next_b(double b0, double r) { return b0 + r; } + + /// Returns the ratio between bin widths. + double r() const noexcept { return r_; } + +private: + // Private constructor below is called by the base class. + friend class piece_base>; + + piece_add(int N, double x0, double b0, double r) + : piece_b0>(N, x0, b0), r_(r) {} + + double r_; +}; + +/** Variant piece + +*/ +template +class piece_variant { + using value_type = Value; + using internal_value_type = detail::get_scale_type; + using piece_variant_type = + boost::variant2::variant, piece_uniform, + piece_multiply, piece_add, + piece_variable>; + +public: + template + explicit piece_variant(const T& p) : piece_(p) {} + + internal_value_type forward(internal_value_type x) const noexcept { + return boost::variant2::visit([x](const auto& p) { return p.forward(x); }, piece_); + } + + internal_value_type inverse(internal_value_type x) const noexcept { + return boost::variant2::visit([x](const auto& p) { return p.inverse(x); }, piece_); + } + + index_type size() const noexcept { + return boost::variant2::visit([](const auto& p) { return p.size(); }, piece_); + } + + internal_value_type x0() const noexcept { + return boost::variant2::visit([](const auto& p) { return p.x0(); }, piece_); + } + + internal_value_type xN() const noexcept { + return boost::variant2::visit([](const auto& p) { return p.xN(); }, piece_); + } + + internal_value_type calc_bin_width(double i) const noexcept { + return boost::variant2::visit([i](const auto& p) { return p.calc_bin_width(i); }, + piece_); + } + + internal_value_type calc_bin_width_last() const noexcept { + return boost::variant2::visit([](const auto& p) { return p.calc_bin_width_last(); }, + piece_); + } + + template + bool has_variant() const noexcept { + return boost::variant2::holds_alternative(piece_); + } + + const piece_variant_type& operator()() const noexcept { return piece_; } + +private: + piece_variant_type piece_; +}; + +} // namespace axis +} // namespace histogram +} // namespace boost + +#endif diff --git a/include/boost/histogram/axis/piecewise.hpp b/include/boost/histogram/axis/piecewise.hpp index aedd12ee..a131154a 100644 --- a/include/boost/histogram/axis/piecewise.hpp +++ b/include/boost/histogram/axis/piecewise.hpp @@ -7,17 +7,22 @@ #ifndef BOOST_HISTOGRAM_AXIS_PIECEWISE_HPP #define BOOST_HISTOGRAM_AXIS_PIECEWISE_HPP +#include #include #include #include #include #include +#include +#include #include +#include #include #include #include #include #include +#include #include #include #include @@ -26,754 +31,219 @@ #include #include -#include - namespace boost { namespace histogram { - namespace axis { -/// Class for starting point -class x_start { -public: - x_start(double value) : value_{value} { - if (!std::isfinite(value)) - BOOST_THROW_EXCEPTION(std::invalid_argument("x_start must be finite")); - } +/** - double value() const { return value_; } +Solution overview: -private: - double value_; -}; +This 1D piecewise axis implementation has two layers. The first layer has to +do with creating immutable transformations from input space X to bin space Y. +The second layer modifies the behavior of these transformation in order to make +them behave like an axis. If axis growth is allowed, this second layer is mutable. +I explain how this works on a high level. -/// Class for stopping point -class x_stop { -public: - x_stop(double value) : value_{value} { - if (!std::isfinite(value)) - BOOST_THROW_EXCEPTION(std::invalid_argument("x_stop must be finite")); - } +Layer 1: Transformation +Transformation pieces: - double value() const { return value_; } +Suppose we want to divide the input space X into 4 bins of unit spacing starting +at x=3. -private: - double value_; -}; +bin 0 1 2 3 + |---|---|---|---| + x=3 x=7 -/// Class for width_start -class width_start { -public: - width_start(double value) : value_{value} { - if (!(0 < value)) - BOOST_THROW_EXCEPTION(std::invalid_argument("width_start must be > 0")); - } +The transformation is described by three pieces of information + forward transformation: y = f(x) = x - 3 + inverse transformation: x = f⁻¹(y) = y + 3 + number of bins: N = 4 - double value() const { return value_; } +There are five supported transformation types: + 1. unit (bin spacing is 1) + 2. uniform (bin spacing is constant) + 3. multiply (bin spacing is multiplied by a constant) + 4. add (bin spacing is added to by a constant) + 5. arbitrary (bin bounds are user supplied points) -private: - double value_; -}; -/** Class that stores the bin spacing option and constant value +Combining transformations: - There are three options for varying bin sizes. +Suppose we wanted to combine the following two transformations. - Uniform |------|------|------|------| - (bi+1 = bi = constant) bi bi+1 + bin 0 1 2 3 bin 0 1 2 3 4 + |---|---|---|---| |-|--|---|----|-----| + x=3 x=7 x=9 x=15 + Unit spacing Add spacing - Multiply |---|-------|---------------| - (bi+1 = bi * constant) bi bi+1 +Typically one wouldn't want to combine two transformations with a gap between them, +but this possibility is supported by the implementation. A piecewise transformation +stores each of these unmodified transformations internally. The external behavior +is that of transformation with bins 0 to 8 shown below. - Add |---|-----|-------|---------| - (bi+1 = bi + constant) bi bi+1 -*/ -class bin_trans { -public: - enum class enum_space { e_uniform, e_multiply, e_add }; + bin 0 1 2 3 nan 4 5 6 7 8 + |---|---|---|---| |-|--|---|----|-----| + x=3 x=7 x=9 x=15 - static bin_trans uniform(double value) { - return bin_trans(enum_space::e_uniform, value); - } - static bin_trans multiply(double value) { - return bin_trans(enum_space::e_multiply, value); - } - static bin_trans add(double value) { return bin_trans(enum_space::e_add, value); } - bin_trans mirror() const noexcept { - if (is_uniform()) return uniform(value_); - if (is_multiply()) return multiply(1 / value_); - if (is_add()) return add(-value_); - assert(false); // TODO: make unnecessary - } +Layer 2: Axis behavior +Bin shifting: - bool is_uniform() const { return spacing_ == enum_space::e_uniform; } - bool is_multiply() const { return spacing_ == enum_space::e_multiply; } - bool is_add() const { return spacing_ == enum_space::e_add; } +Suppose we wanted to use the transformation above in a growable axis. This axis starts +off with 9 bins. Two kinds of bins can be added: underflow and overflow. Underflow bins +are added below x=3. Overflow bins are added above x=15. Adding bins between x=7 and +x=9 is not supported. The mutable transformation `x_bin_shift` accumulates: + - added underflow bins, and + - added overflow bins. +Using this information, x_bin_shift augments the behavior of the contained piecewise +transformation as follows. The size (i.e., number of bins) of the transformation is +increased by the number of underflow and overflow bins. The bin values (i.e., Y axis) +is shifted by the number of underflow bins in the forward and inverse calculations. - enum_space spacing() const { return spacing_; } - double value() const { return value_; } +Linear axis: -private: - bin_trans(enum_space spacing, double value) : spacing_{spacing}, value_{value} {} +Wrapping the, possibly shifted, transformation in an `int_resolver_linear` makes the +transformation behave like an axis. This linear int resolver implemenets an index method. - enum_space spacing_{enum_space::e_uniform}; - double value_; -}; +Circular axis: -/** Class for a piece of a piecewise axis +Suppose we wanted to make the above described piecewise transformation circular with +a relevant period x=3 to x=23. This transformaion is copied below for convenience. - Each piece of a piecewise axis has a start point x0 and an end point xN. - It is required that x0 < xN. + bin 0 1 2 3 nan 4 5 6 7 8 + |---|---|---|---| |-|--|---|----|-----| + x=3 x=7 x=9 x=15 - |--------------------------------------| - x0 xN +To do this, we wrap the transformation in an `int_resolver_circular` with the bounds +x=3 and x=23. This circular transform augments the behavior of the contained +transformation. The input (i.e., the X axis) is wrapped to the domain x=3 to x=23 +before being given to the contained transformation. This circular int resolver +implemenets an index method allowing it to act as an axis. - Each piece is divided into n bins. The bin spacing is determined by the bin spacing - option and value. +*/ - There are two ways to create a piece: rollout and force. With the rollout method, bins - are added one at a time. If the bin number is positive (n > 0), n bins are added to the - right. If the bin number is negative (n < 0), |n| bins are added to the left. With the - force method, one gives the interval x0 to xN and the number of bins n. The problem - finds the bin width. +/** Piecewise -- putting pieces together */ -template -class piece { +template +class piecewise { using value_type = Value; - using unit_type = detail::get_unit_type; using internal_value_type = detail::get_scale_type; public: - /** Creates a piece with the rollout method - - @param n number of bins (right if positive, left if negative) - @param x_start starting point - @param width_start starting bin width - @param bt bin spacing option and value - - The constructor throws `std::invalid_argument` if the bin number is zero or if the - bin spacing option is invalid. - - Example: Start at x = 1 with an initial bin size of 1. Rollout 4 bins to the right, - doubling the bin size each time. - |-|---|-------|---------------| - 0 2 4 6 8 10 12 14 16 - - */ - static piece rollout(int n, x_start x, width_start b0, bin_trans bt) { - assert(n != 0); - if (0 < n) { - return rollout_right(n, x, b0, bt); - } else { - return rollout_left(n, x, b0, bt); - } - } - - /** Creates a piece with the force method - - @param n number of bins - @param x_start starting point - @param x_stop stopping point - @param bt bin spacing option and value - - The constructor throws `std::invalid_argument` if: - 1. the bin number is zero, - 2. x_stop <= x_start, or - 3. the bin spacing option is invalid. - - Example: Input interval [1, 16] - - |-----------------------------| - 0 2 4 6 8 10 12 14 16 - - Place 4 bins. Each bin is double the size of the previous bin. - - |-|---|-------|---------------| - 0 2 4 6 8 10 12 14 16 - */ - static piece force(int n, x_start x_start, x_stop x_stop, bin_trans bt) { - // TODO: Check that x_start < x_stop - const bool is_ordered = x_start.value() < x_stop.value(); - if (!is_ordered) - BOOST_THROW_EXCEPTION(std::invalid_argument("x_start must be < x_stop")); - - const width_start b0_obj = piece::calc_b0(x_start, x_stop, n, bt); - return piece(n, x_start, b0_obj, bt); + /// Creates a piecewise starting with an initial piece. + template + explicit piecewise(const T& p) { + v_pieces_.push_back(PieceType(p)); } + /// The mapping from input space X to bin space Y template T forward(T x) const noexcept { - static_assert(std::is_floating_point::value, "T must be a floating point type"); - // Runs in hot loop, please measure impact of changes - T f_x; - if (bt_.is_uniform()) { - f_x = forward_uniform(x); - } else if (bt_.is_multiply()) { - f_x = forward_multiply(x); - } else if (bt_.is_add()) { - f_x = forward_add(x); - } else { - assert(false); // TODO: make unnecessary - } - return f_x + accumulated_left_shift_; - } - index_type size() const noexcept { - return size_ic_ + accumulated_left_shift_ + accumulated_right_shift_; - } - - double x0() const noexcept { return x0_; } - double b0() const noexcept { return b0_; } - bin_trans bt() const noexcept { return bt_; } - index_type size_ic() const noexcept { return size_ic_; } - - internal_value_type xN() const noexcept { return xN_; } + // First piece + const auto& piece_0 = v_pieces_.front(); + if (x < piece_0.xN()) return piece_0.forward(x); - index_type accumulated_left_shift() const noexcept { return accumulated_left_shift_; } - index_type accumulated_right_shift() const noexcept { return accumulated_right_shift_; } + // Middle pieces + index_type offset = piece_0.size(); + const int n_middle = v_pieces_.size() - 1; + for (int j = 1; j < n_middle; ++j) { + const auto& p = v_pieces_[j]; - // Bin ending at x_i has width - double calc_bin_width(double i) const noexcept { - const auto x_i = reverse(i); - const auto x_i_prev = reverse(i - 1); - return x_i - x_i_prev; - } + if (x < p.x0()) { + return std::numeric_limits::quiet_NaN(); // Between discontinuous pieces + } else if (x < p.xN()) { + return offset + p.forward(x); + } - static width_start calc_b0(x_start x0, x_stop xN, int N, bin_trans bt) { - double b0; - if (bt.is_uniform()) { - b0 = calc_b0_uniform(x0, xN, N); - } else if (bt.is_multiply()) { - b0 = calc_b0_multiply(x0, xN, N, bt); - } else if (bt.is_add()) { - b0 = calc_b0_add(x0, xN, N, bt); - } else { - assert(false); // TODO: make unnecessary + offset += p.size(); } - return width_start(b0); - } - template - internal_value_type reverse(T x) const noexcept { - if (bt_.is_uniform()) { - return reverse_uniform(x); - } else if (bt_.is_multiply()) { - return reverse_multiply(x); - } else if (bt_.is_add()) { - return reverse_add(x); + const auto& piece_N = get_piece_R(); + if (piece_N.x0() <= x) { + return offset + get_piece_R().forward(x); // Last piece } else { - assert(false); // TODO: make unnecessary + return std::numeric_limits::quiet_NaN(); // Before last piece } } -private: - static piece rollout_right(int n, x_start xL, width_start b0, bin_trans bt) { - assert(0 < n); - return piece(n, xL, b0, bt); - } - - static piece rollout_left(int n_neg, x_start xR, width_start b0_R, - bin_trans bt_right) { - assert(n_neg < 0); - const int n_pos = -n_neg; - - // Roll out right - const piece piece_right(n_pos, xR, b0_R, bt_right); - - // Get width of last bin of right piece - const double width_last_bin = piece_right.calc_bin_width(n_pos); + /** Extrapolates the bin spacing of the right side of the piecewise axis. - // Get location of last bin of right piece - const double x_last_bin = piece_right.reverse(n_pos); + @param N number of bins to extrapolate + @param args arguments to pass to the piece constructor - // Reflect last bin of right piece to get xL - const double delta = x_last_bin - xR.value(); - const double xL = xR.value() - delta; + Example: you have a piecewise axis with unit spacing and you want to double the + spacing of the right side for 3 bins. The following code will do this. - return piece(n_pos, x_start(xL), width_start(width_last_bin), bt_right.mirror()); - } + this->extrapolate_R(3, 2.0); - // Uniform spacing - // - // b0 b0 b0 b0 - // N=0 N=1 N=2 N=3 N=4 - // |------|------|------|------| - // x0 x4 - // - // The sequence - // x0 = x0 - // x1 = x0 + b0 - // x2 = x0 + b0 * 2 - // x3 = x0 + b0 * 3 - // ... - // xN = x0 + b0 * N - // - // Formula for N in terms of xN - // N = (xN - x0) / b0 - // - template - internal_value_type forward_uniform(T x) const noexcept { - // Runs in hot loop, please measure impact of changes - return ((x / unit_type{}) - x0_) / b0_; - } - - // Starting with: - // xN = x0 + b0 * N - // Solve for b0 to get: - // b0 = (xN - x0) / N - static double calc_b0_uniform(x_start x0, x_stop xN, int N) { - return (xN.value() - x0.value()) / N; - } - - template - internal_value_type reverse_uniform(T x) const noexcept { - return x0_ + x * b0_; - } - - // Multiply spacing - // - // b0 b0*r b0*r² b0*r³ - // N=0 N=1 N=2 N=3 N=4 - // |---|-----|-----------|------------------| - // x0 x4 - // - // The sequence (for some value of ψ) - // x0 = 1 * b0 / (r - 1) + ψ - // x1 = r * b0 / (r - 1) + ψ - // x2 = r² * b0 / (r - 1) + ψ - // x3 = r³ * b0 / (r - 1) + ψ - // ... - // xN = rᴺ * b0 / (r - 1) + ψ - // - // Note: the first bin spacing is b0 - // x1 - x0 = r * b0 / (r - 1) - 1 * b0 / (r - 1) - // = (r - 1) * b0 / (r - 1) - // = b0 - // - // Find a formula for N - // xN - x0 = b0 / (r - 1) * rᴺ - b0 / (r - 1) * 1 - // = b0 / (r - 1) * (rᴺ - 1) - // - // N in terms of xN - // N = log2(((xN - x0) / b0) * (r - 1) + 1) / log2(r) - // - template - internal_value_type forward_multiply(T x) const noexcept { - // Runs in hot loop, please measure impact of changes - const auto z = forward_uniform(x); // z = (x - x0) / b0 - const auto r = bt_.value(); - return std::log2(z * (r - 1) + 1) / std::log2(r); - } - - // Starting with: - // xN - x0 = b0 / (r - 1) * (rᴺ - 1) - // Solve for b0 to get: - // b0 = (xN - x0) * (r - 1) / (rᴺ - 1) - static double calc_b0_multiply(x_start x0, x_stop xN, int N, bin_trans bt) { - const auto r = bt.value(); - return (xN.value() - x0.value()) * (r - 1) / (std::pow(r, N) - 1); - } - - template - internal_value_type reverse_multiply(T x) const noexcept { - const auto r = bt_.value(); - return x0_ + b0_ * (std::pow(r, x) - 1) / (r - 1); + <- - - -new piece- - - - -> + |-|-|-|-|-|---|-------|---------------| + 3 4 5 6 7 8 9 10 14 22 + */ + template + void extrapolate_R(int N, Args... args) { + double x0 = xN(); // Right side of the right piece becomes x0 for the new piece + double b_prev = get_piece_R().calc_bin_width_last(); + double b0 = P::next_b(b_prev, args...); + const auto p = P::create(N, x0, b0, args...); + add_R(p); } - // Add spacing - // - // b0 b0+r b0+2r b0+3r - // N=0 N=1 N=2 N=3 N=4 - // |---|-----|-----------|------------------| - // x0 x4 - // - // The sequence - // x0 = x0 - // x1 = x0 + 1 * b0 - // x2 = x0 + 2 * b0 + r - // x3 = x0 + 3 * b0 + r * (1 + 2) - // x4 = x0 + 4 * b0 + r * (1 + 2 + 3) - // ... - // x = x0 + N * b0 + r * (1 + 2 + ... + N-1) - // = x0 + N * b0 + r * (N * (N - 1) / 2) - // - // Multiply by 2 - // 2 * x = 2 * x0 + 2 * N * b0 + r * (N * (N - 1)) - // Expand terms: - // 2 * x = 2 * x0 + 2 * N * b0 + rN² - rN - // Put all terms on one side: - // 0 = 2 * x0 + 2 * N * b0 + rN² - rN - 2 * x - // Solve quadratic equation and take positive root - // N = (-1 + sqrt(1 + 8 * (x - x0) / b0 * r)) / 2 - // - template - internal_value_type forward_add(T x) const noexcept { - // Runs in hot loop, please measure impact of changes - const auto z = forward_uniform(x); // z = (x - x0) / b0 - const auto r = bt_.value(); - return (-1 + std::sqrt(1 + 8 * z * r)) / 2; - } + /** Adds a piece to the right side of the piecewise axis. - // Starting with: - // 2 * x = 2 * x0 + 2 * N * b0 + r * (N * (N - 1)) - // Solve for b0 to get: - // b0 = 2 * (xN - x0) / (N * (N - 1)) - r - // - static double calc_b0_add(x_start x0, x_stop xN, int N, bin_trans bt) { - const auto r = bt.value(); - return 2 * (xN.value() - x0.value()) / (N * (N - 1)) - r; - } + @param p piece to add - template - internal_value_type reverse_add(T x) const noexcept { - const auto r = bt_.value(); - return x0_ + x * b0_ + r * (x * (x - 1)) / 2; - } - - // Private constructor - piece(int n, x_start x0, width_start b0, bin_trans bt) - : size_ic_(n), x0_(x0.value()), b0_(b0.value()), bt_(bt) { - xN_ = reverse(size_ic_); + The constructor throws `std::invalid_argument` if the left side of the new piece + does not exactly match the right side of the right piece. + */ + template + void add_R(const P& p) { + if (p.x0() != xN()) // Right side of the right piece is not equal to the left side of + // the new piece + BOOST_THROW_EXCEPTION(std::invalid_argument("x0 != xN")); + add_R_gap_okay(p); } - index_type size_ic_{-9999}; - internal_value_type x0_{-9999.0}; - double b0_{-9999.0}; - bin_trans bt_{}; - - internal_value_type xN_{-9999.0}; - - index_type accumulated_left_shift_{0}; - index_type accumulated_right_shift_{0}; -}; - -/***Extrapolation and Attachment*** - - Pieces can be added to the right or left side of an existing piece. Pieces are - added with one of two approaches: extrapolation or attachment. Examples of each are - shown below. - Extrapolate 2 bins to the right side, doubling the bin size each time. - 0 2 4 6 8 - |---|---|---|---| - <- - - New piece - - -> - |---|---|---|---|-------|---------------| - 0 2 4 6 8 10 12 14 16 18 20 - Attach a new piece with a uniform bin spacing to the right side. - 0 2 4 6 8 - |---|---|---|---| - <- - - New piece - - -> - |---|---|---|---|-|-|-|-|-|-|-|-|-|-|-|-| - 0 2 4 6 8 10 12 14 16 18 20 - -***Syntax*** - -using P = piecewise; - -auto pa = P( // Create a new piece - P::x_start(1), // Start at x_start = 1 - P::width_start(1), // Start with a bin size of width_start = 1 - P::bin_trans::uniform(1), // Add bins of constant size 1 - P::stop::n_bins(4)); // Add 4 bins - -pa.add_right( // Add to the right side - P::bin_trans::multiply(2), // Extrapolate, doubling the bin size each time - P::stop::x(16)); // Stop when the piece contains x = 16 - -pa.add_left( // Add to the left side - P::width_start(1), // Start with a bin size of width_start = 1 - P::bin_trans::add(1), // Add bins of constant size 1 - P::stop::x(4)); // Stop when the piece contains x = 4 -*/ -template -class piecewise { - using value_type = Value; - using unit_type = detail::get_unit_type; - using internal_value_type = detail::get_scale_type; - -public: - /// Construct a piecewise axis with an initial piece. - explicit piecewise(piece p) { v_pieces_.push_back(p); } + /** Adds a piece to the right side of the piecewise axis. - /// - template - T forward(T x) const noexcept { - // Runs in hot loop, please measure impact of changes - index_type offset = 0; - const int n_pieces = v_pieces_.size(); - for (int j = 0; j < n_pieces; ++j) { - const auto& p = v_pieces_[j]; - if (x < p.x1 || j == n_pieces - 1) return offset + p.index(x); - } - } - - /// Shifts the axis - void shift_axis(index_type n) { - if (n < 0) { - v_pieces_.front().shift_axis(n); - } else if (0 < n) { - // Get last element of v_pieces - v_pieces_.back().shift_axis(n); - } - } + @param p piece to add - void add_right(int n, bin_trans bt) { - double width_start = get_right_width_start(); - add_right(n, width_start, bt); - } - void add_right(int n, width_start width_start, bin_trans bt) { - double x_start = get_right_x(); + The constructor throws `std::invalid_argument` if the left side of the new piece + is less than the right side of the right piece. + */ + template + void add_R_gap_okay(const P& p) { + // Check that the left side of the new piece is equal to or greater than the right + // side of the right piece + if (p.x0() < xN()) BOOST_THROW_EXCEPTION(std::invalid_argument("x0 < xN")); - const auto p = piece::rollout(n, x_start, width_start, bt); - v_pieces_.push_back(p); + v_pieces_.push_back(PieceType(p)); } + /// The number of bins in the piecewise. index_type size() const noexcept { index_type size = 0; for (const auto& p : v_pieces_) { size += p.size(); } return size; - // One liner for the above function: - // return std::accumulate(v_pieces_.begin(), v_pieces_.end(), 0); - } - -private: - const piece& get_left() const noexcept { - assert(!v_pieces_.empty()); - return v_pieces_.front(); - } - - const piece& get_right() const noexcept { - assert(!v_pieces_.empty()); - return v_pieces_.back(); } - double get_left_x() const noexcept { return get_left().x0; } - double get_right_x() const noexcept { return get_right().x1; } + /// The location of the right side of the right piece. + double xN() const noexcept { return get_piece_R().xN(); } - double get_left_width_start() const noexcept { return get_left().width_start; } - double get_right_width_start() const noexcept { return get_right().width_start; } - - std::vector> v_pieces_{}; -}; - -/// -template -class int_shifter { - using value_type = Value; - -public: - /// Shifts the axis - void shift_axis(index_type n) { - if (n < 0) { - // Need to offset by this in the future - accumulated_left_shift_ += std::abs(n); - } else if (0 < n) { - // Need to offset by this in the future - accumulated_right_shift_ += n; - } - } + /// The location of the left side of the left piece. + double x0() const noexcept { return v_pieces_.front().x0(); } private: - Payload payload_; - index_type accumulated_left_shift_{0}; - index_type accumulated_right_shift_{0}; -}; - -/// -template -class int_resolver { - using value_type = Value; - -public: - index_type index(value_type x) const noexcept { - // Runs in hot loop, please measure impact of changes - - y = payload_.index(x); - - if (y < size()) { - if (0 <= y) - return static_cast(y); // 0 <= i < size - else - return -1; // i < 0 - } - - // upper edge of last bin is inclusive if overflow bin is not present - if constexpr (std::is_floating_point::value) { - if (!options_type::test(option::overflow) && y == size()) return size() - 1; - } - - return size(); // also returned if x is NaN - } - - std::pair update(value_type x) noexcept { - // Runs in hot loop, please measure impact of changes - - y = payload_.index(x); - - if (y < size()) { - if (0 <= y) { - const auto i_int = static_cast(y); - return {i_int, 0}; - } else if (y != -std::numeric_limits::infinity()) { - const auto i_int = static_cast(std::floor(y)); - shift_axis(i_int); - return {0, -i_int}; - } else { - // i is -infinity - return {-1, 0}; - } - } - // i either beyond range, infinite, or NaN - if (y < std::numeric_limits::infinity()) { - const auto i_int = static_cast(y); - const auto n = i_int - size() + 1; - shift_axis(n); - return {y, -n}; - } - // z either infinite or NaN - return {size(), 0}; - } - - index_type size() const noexcept { return payload_.size(); } - -private: - Payload payload_; -}; - -/// -template -class int_resolver_circular { - using value_type = Value; - -public: - index_type index(value_type x) const noexcept { - // Runs in hot loop, please measure impact of changes - - y = payload_.index(x); - - // Not finite --> overflow bin - if constexpr (std::is_floating_point::value) { - if (!std::isfinite(y)) return payload_.size(); - } - - return y % payload_.size(); - } - -private: - Payload payload_; -}; - -/** Axis for unit intervals on the real line. - - The most common binning strategy. Very fast. Binning is a O(1) operation. - - If the axis has an overflow bin (the default), a value on the upper edge of the last - bin is put in the overflow bin. The axis range represents a semi-open interval. - - If the overflow bin is deactivated, then a value on the upper edge of the last bin is - still counted towards the last bin. The axis range represents a closed interval. - - The options `growth` and `circular` are mutually exclusive. - - @tparam Value input value type, must be floating point. - @tparam Transform builtin or user-defined transform type. - @tparam MetaData type to store meta data. - @tparam Options see boost::histogram::axis::option. - */ -template -class unit_regular : public iterator_mixin>, - // protected detail::replace_default, - public metadata_base_t { - // these must be private, so that they are not automatically inherited - using value_type = Value; - using metadata_base = metadata_base_t; - using metadata_type = typename metadata_base::metadata_type; - using options_type = - detail::replace_default; - - using unit_type = detail::get_unit_type; - using internal_value_type = detail::get_scale_type; - using piecewise_type = piecewise; - -public: - constexpr unit_regular() = default; - - unit_regular(piecewise pw, metadata_type meta = {}, - options_type options = {}) - : metadata_base(std::move(meta)), piecewise_(std::move(pw)) { - // static_asserts were moved here from class scope to satisfy deduction in gcc>=11 - static_assert(std::is_nothrow_move_constructible::value, - "piecewise must be no-throw move constructible"); - static_assert(std::is_nothrow_move_assignable::value, - "piecewise must be no-throw move assignable"); - static_assert(std::is_floating_point::value, - "unit_regular axis requires floating point type"); - static_assert(!(options.test(option::circular) && options.test(option::growth)), - "circular and growth options are mutually exclusive"); - if (size() <= 0) BOOST_THROW_EXCEPTION(std::invalid_argument("bins > 0 required")); - } - - /// Constructor used by algorithm::reduce to shrink and rebin (not for users). - unit_regular(const unit_regular& src, index_type begin, index_type end, unsigned merge) - : unit_regular(src.transform(), (end - begin) / merge, src.value(begin), - src.value(end), src.metadata()) { - assert(false); - assert((end - begin) % merge == 0); - if (options_type::test(option::circular) && !(begin == 0 && end == src.size())) - BOOST_THROW_EXCEPTION(std::invalid_argument("cannot shrink circular axis")); - } - - /// Return index for value argument. - index_type index(value_type x) const noexcept { return piecewise_.index(x); } - - /// Returns index and shift (if axis has grown) for the passed argument. - std::pair update(value_type x) noexcept { - assert(options_type::test(option::growth)); - return piecewise_.update(x); - } - - /// Return value for fractional index argument. - value_type value(real_index_type i) const noexcept { - assert(false); - return static_cast(i); - } - - /// Return bin for index argument. - decltype(auto) bin(index_type idx) const noexcept { - return interval_view(*this, idx); + // Gets the right-most piece. + const PieceType& get_piece_R() const noexcept { + assert(!v_pieces_.empty()); + return v_pieces_.back(); } - /// Returns the number of bins, without over- or underflow. - index_type size() const noexcept { return size_; } - - /// Returns the options. - static constexpr unsigned options() noexcept { return options_type::value; } - - // template - // bool operator==(const unit_regular& o) const noexcept { - // return detail::relaxed_equal{}(transform(), o.transform()) && size() == o.size() && - // min_ == o.min_ && delta_ == o.delta_ && - // detail::relaxed_equal{}(this->metadata(), o.metadata()); - // } - // template - // bool operator!=(const unit_regular& o) const noexcept { - // return !operator==(o); - // } - - // template - // void serialize(Archive& ar, unsigned /* version */) { - // ar& make_nvp("transform", static_cast(*this)); - // ar& make_nvp("size", size_); - // ar& make_nvp("meta", this->metadata()); - // ar& make_nvp("min", min_); - // ar& make_nvp("delta", delta_); - // } - -private: - index_type size_{0}; - internal_value_type min_{0}; - piecewise piecewise_{}; + std::vector v_pieces_{}; }; } // namespace axis diff --git a/test/CMakeLists.txt b/test/CMakeLists.txt index cc4b0734..f208a82a 100644 --- a/test/CMakeLists.txt +++ b/test/CMakeLists.txt @@ -54,7 +54,9 @@ boost_test(TYPE run SOURCES axis_category_test.cpp) boost_test(TYPE run SOURCES axis_integer_test.cpp) boost_test(TYPE run SOURCES axis_option_test.cpp) boost_test(TYPE run SOURCES axis_regular_test.cpp) +boost_test(TYPE run SOURCES axis_piece_test.cpp) boost_test(TYPE run SOURCES axis_piecewise_test.cpp) +boost_test(TYPE run SOURCES axis_int_resolver_test.cpp) boost_test(TYPE run SOURCES axis_traits_test.cpp) boost_test(TYPE run SOURCES axis_variable_test.cpp) boost_test(TYPE run SOURCES axis_variant_test.cpp) diff --git a/test/Jamfile b/test/Jamfile index cef1a81c..7523d39e 100644 --- a/test/Jamfile +++ b/test/Jamfile @@ -57,7 +57,9 @@ alias cxx14 : [ run axis_integer_test.cpp ] [ run axis_option_test.cpp ] [ run axis_regular_test.cpp ] + [ run axis_piece_test.cpp ] [ run axis_piecewise_test.cpp ] + [ run axis_int_resolver_test.cpp ] [ run axis_traits_test.cpp ] [ run axis_variable_test.cpp ] [ run axis_variant_test.cpp ] diff --git a/test/axis_int_resolver_test.cpp b/test/axis_int_resolver_test.cpp new file mode 100644 index 00000000..d40b19c4 --- /dev/null +++ b/test/axis_int_resolver_test.cpp @@ -0,0 +1,111 @@ +// CopyR 2015-2019 Hans Dembinski +// +// Distributed under the Boost Software License, Version 1.0. +// (See accompanying file LICENSE_1_0.txt +// or copy at http://www.boost.org/LICENSE_1_0.txt) + +#include +#include +#include +#include +#include +#include +#include +#include +#include "axis.hpp" +#include "is_close.hpp" +#include "std_ostream.hpp" +#include "str.hpp" +#include "throw_exception.hpp" + +#include +#include +#include + +using namespace boost::histogram; + +int main() { + // int_resolver_circular + { + using type_multiply = axis::piece_multiply; + + // Has x0 = 4, x1 = 5, x2 = 7, x3 = 11, x4 = 17 + const int k_N = 3; + const double k_x0 = 3.0; + const double k_b0 = 1.0; + const double k_r = 2.0; + + const auto p = type_multiply::create(k_N, k_x0, k_b0, k_r); + const auto a = axis::int_resolver_circular(p, 2, 22); + + BOOST_TEST_EQ(a.index(3), 0); + BOOST_TEST_EQ(a.index(4), 1); + BOOST_TEST_EQ(a.index(6), 2); + BOOST_TEST_EQ(a.index(10), 3); + BOOST_TEST_EQ(a.index(18), 4); + // BOOST_TEST_EQ(a.index(5), k_x0 + k_b0); + } + + // int_resolver_linear + { + using type_uniform = axis::piece_uniform; + const type_uniform p_uniform = type_uniform::solve_b0(4, -2, +2); + const auto a = axis::int_resolver_linear(p_uniform); + + // Copied from axis_regular_test.cpp + BOOST_TEST_EQ(a.value(0), -2); + BOOST_TEST_EQ(a.value(1), -1); + BOOST_TEST_EQ(a.value(2), 0); + BOOST_TEST_EQ(a.value(3), 1); + BOOST_TEST_EQ(a.value(4), 2); + // BOOST_TEST_EQ(a.bin(-1).lower(), -std::numeric_limits::infinity()); + // BOOST_TEST_EQ(a.bin(-1).upper(), -2); + // BOOST_TEST_EQ(a.bin(a.size()).lower(), 2); + // BOOST_TEST_EQ(a.bin(a.size()).upper(), std::numeric_limits::infinity()); + BOOST_TEST_EQ(a.index(-10.), -1); + BOOST_TEST_EQ(a.index(-2.1), -1); + BOOST_TEST_EQ(a.index(-2.0), 0); + BOOST_TEST_EQ(a.index(-1.1), 0); + BOOST_TEST_EQ(a.index(0.0), 2); + BOOST_TEST_EQ(a.index(0.9), 2); + BOOST_TEST_EQ(a.index(1.0), 3); + BOOST_TEST_EQ(a.index(10.), 4); + BOOST_TEST_EQ(a.index(-std::numeric_limits::infinity()), -1); + BOOST_TEST_EQ(a.index(std::numeric_limits::infinity()), 4); + BOOST_TEST_EQ(a.index(std::numeric_limits::quiet_NaN()), 4); + } + + // with growth + { + using pii_t = std::pair; + using type_uniform = axis::piece_uniform; + using type_shift = axis::bin_shift_integrator; + using type_resolver = axis::int_resolver_linear; + + const auto p_uniform = type_uniform::solve_b0(1, 0, 1); + const auto p_shift = type_shift(p_uniform); + auto a = type_resolver(p_shift); + + // Copied from axis_regular_test.cpp + BOOST_TEST_EQ(a.size(), 1); + BOOST_TEST_EQ(a.update(0), pii_t(0, 0)); + BOOST_TEST_EQ(a.size(), 1); + BOOST_TEST_EQ(a.update(1), pii_t(1, -1)); + BOOST_TEST_EQ(a.size(), 2); + BOOST_TEST_EQ(a.value(0), 0); + BOOST_TEST_EQ(a.value(2), 2); + BOOST_TEST_EQ(a.update(-1), pii_t(0, 1)); + BOOST_TEST_EQ(a.size(), 3); + BOOST_TEST_EQ(a.value(0), -1); + BOOST_TEST_EQ(a.value(3), 2); + BOOST_TEST_EQ(a.update(-10), pii_t(0, 9)); + BOOST_TEST_EQ(a.size(), 12); + BOOST_TEST_EQ(a.value(0), -10); + BOOST_TEST_EQ(a.value(12), 2); + BOOST_TEST_EQ(a.update(std::numeric_limits::infinity()), pii_t(a.size(), 0)); + BOOST_TEST_EQ(a.update(std::numeric_limits::quiet_NaN()), pii_t(a.size(), 0)); + BOOST_TEST_EQ(a.update(-std::numeric_limits::infinity()), pii_t(-1, 0)); + } + + return boost::report_errors(); +} diff --git a/test/axis_piece_test.cpp b/test/axis_piece_test.cpp new file mode 100644 index 00000000..0fecb153 --- /dev/null +++ b/test/axis_piece_test.cpp @@ -0,0 +1,213 @@ +// CopyR 2015-2019 Hans Dembinski +// +// Distributed under the Boost Software License, Version 1.0. +// (See accompanying file LICENSE_1_0.txt +// or copy at http://www.boost.org/LICENSE_1_0.txt) + +#include +#include +#include +#include +#include +#include +#include "axis.hpp" +#include "is_close.hpp" +#include "std_ostream.hpp" +#include "str.hpp" +#include "throw_exception.hpp" + +#include +#include +#include + +using namespace boost::histogram; + +// +class tester_piece_b0 { +public: + static tester_piece_b0 case_1() { return tester_piece_b0(4, -1.1, 1.2); } + + tester_piece_b0(int bins, double x0, double b0) : size_(bins), x0_(x0), b0_(b0) {} + + template + PieceType create(Args... args) const { + return PieceType::create(size_, x0_, b0_, args...); + } + + template + void test_create(const PieceType& p, Lambda&& lambda) const { + + // Check abstract piece member data + BOOST_TEST_EQ(p.size(), size_); + BOOST_TEST_EQ(p.x0(), x0_); + // xN is tested below + + double bin_width = b0_; + double x = x0_; + for (int i = 0; i < size_; ++i) { + // Check that forward forward transformation (x --> float_index) is correct + BOOST_TEST_IS_CLOSE(p.forward(x), i, 1.0e-12); + x += bin_width; + + // Check that the bin width is correct + BOOST_TEST_IS_CLOSE(p.calc_bin_width(i), bin_width, 1.0e-12); + bin_width = lambda(bin_width); + // tester_next_width(p, &bin_width, args...); + } + + // Check xN + BOOST_TEST_IS_CLOSE(p.xN(), x, 1.0e-12); + + // Check that inverse transformation is correct + for (double v = -9.0; v < 9.0; v += 0.01) { + BOOST_TEST_IS_CLOSE(p.inverse(p.forward(v)), v, 1.0e-12); + } + } + + int size() const { return size_; } + double x0() const { return x0_; } + double b0() const { return b0_; } + +private: + int size_; + double x0_; + double b0_; +}; + +template +void test_equal(const PieceType& p1, const PieceType& p2) { + BOOST_TEST_EQ(p1.size(), p2.size()); + BOOST_TEST_IS_CLOSE(p1.x0(), p2.x0(), 1.0e-12); + BOOST_TEST_IS_CLOSE(p1.xN(), p2.xN(), 1.0e-12); + + // Checks (indirectly) that the parameters of the piece are the same + for (int i = 0; i < p1.size(); ++i) { + BOOST_TEST_IS_CLOSE(p1.calc_bin_width(i), p2.calc_bin_width(i), 1.0e-12); + } +} + +template +void test_solve_b0(const PieceType& p_R, Args... args) { + // Construct force from create + const auto p2 = PieceType::solve_b0(p_R.size(), p_R.x0(), p_R.xN(), args...); + + test_equal(p_R, p2); +} + +int main() { + BOOST_TEST(std::is_nothrow_move_assignable>::value); + BOOST_TEST(std::is_nothrow_move_constructible>::value); + + // Test unit + { + using p_type = axis::piece_unit; + const auto p = p_type::create(4, 2); + + BOOST_TEST_EQ(p.size(), 4); + BOOST_TEST_EQ(p.x0(), 2); + BOOST_TEST_EQ(p.xN(), 6); + + BOOST_TEST_EQ(p.forward(2), 0); + BOOST_TEST_EQ(p.forward(3), 1); + BOOST_TEST_EQ(p.forward(4), 2); + BOOST_TEST_EQ(p.forward(5), 3); + BOOST_TEST_EQ(p.forward(6), 4); + + BOOST_TEST_EQ(p.inverse(0), 2); + BOOST_TEST_EQ(p.inverse(1), 3); + BOOST_TEST_EQ(p.inverse(2), 4); + BOOST_TEST_EQ(p.inverse(3), 5); + BOOST_TEST_EQ(p.inverse(4), 6); + } + + // Test arbitrary + { + using p_type = axis::piece_variable; + const auto p = p_type::create(std::vector{1.1, 2.6, 3.8, 4.2, 5.9}); + + BOOST_TEST_EQ(p.size(), 4); + BOOST_TEST_EQ(p.x0(), 1.1); + BOOST_TEST_EQ(p.xN(), 5.9); + + BOOST_TEST_EQ(p.forward(1.1), 0); + BOOST_TEST_EQ(p.forward(2.6), 1); + BOOST_TEST_EQ(p.forward(3.8), 2); + BOOST_TEST_EQ(p.forward(4.2), 3); + // BOOST_TEST_EQ(p.forward(5.9), 4); // Test broken + + BOOST_TEST_EQ(p.inverse(0), 1.1); + BOOST_TEST_EQ(p.inverse(1), 2.6); + BOOST_TEST_EQ(p.inverse(2), 3.8); + BOOST_TEST_EQ(p.inverse(3), 4.2); + BOOST_TEST_EQ(p.inverse(4), 5.9); + } + + // Test b0 pieces + { + const tester_piece_b0 p1 = tester_piece_b0::case_1(); + + // Piece uniform + { + using p_type = axis::piece_uniform; + const auto p = p1.create(); + auto lambda_b_next = [](double b) { return b; }; + { + p1.test_create(p, lambda_b_next); + test_solve_b0(p); + } + + // Variant + { + const auto p_orig = p1.create(); + const auto p = axis::piece_variant(p_orig); + BOOST_TEST(p.has_variant>()); + BOOST_TEST(!p.has_variant>()); + p1.test_create(p, lambda_b_next); + } + } + + // Piece multiply + { + using p_type = axis::piece_multiply; + const double k_r_multiply = 1.1; + auto lambda_b_next = [&](double b) { return b * k_r_multiply; }; + { + const auto p = p1.create(k_r_multiply); + p1.test_create(p, lambda_b_next); + test_solve_b0(p, k_r_multiply); + } + + // Variant + { + const auto p_orig = p1.create(k_r_multiply); + const auto p = axis::piece_variant(p_orig); + BOOST_TEST(p.has_variant>()); + BOOST_TEST(!p.has_variant>()); + p1.test_create(p, lambda_b_next); + } + }; + + // Piece add + { + using p_type = axis::piece_add; + const double k_r_add = 0.05; + auto lambda_b_next = [&](double b) { return b + k_r_add; }; + { + const auto p = p1.create(k_r_add); + p1.test_create(p, lambda_b_next); + test_solve_b0(p, k_r_add); + } + + // Variant + { + const auto p_orig = p1.create(k_r_add); + const auto p = axis::piece_variant(p_orig); + BOOST_TEST(p.has_variant>()); + BOOST_TEST(!p.has_variant>()); + p1.test_create(p, lambda_b_next); + } + } + } + + return boost::report_errors(); +} diff --git a/test/axis_piecewise_test.cpp b/test/axis_piecewise_test.cpp index c35f4455..bf3d73f6 100644 --- a/test/axis_piecewise_test.cpp +++ b/test/axis_piecewise_test.cpp @@ -1,4 +1,4 @@ -// Copyright 2015-2019 Hans Dembinski +// CopyR 2015-2019 Hans Dembinski // // Distributed under the Boost Software License, Version 1.0. // (See accompanying file LICENSE_1_0.txt @@ -6,6 +6,7 @@ #include #include +#include #include #include #include @@ -20,373 +21,57 @@ #include #include -int main() { - using namespace boost::histogram; - // using def = use_default; - // namespace tr = axis::transform; - // namespace op = axis::option; - - // BOOST_TEST(std::is_nothrow_move_assignable>::value); - // BOOST_TEST(std::is_nothrow_move_constructible>::value); +#include - BOOST_TEST(std::is_nothrow_move_assignable>::value); - BOOST_TEST(std::is_nothrow_move_constructible>::value); +using namespace boost::histogram; - // Construct rollout right +int main() { + // Piecewise Connected { - const auto p = axis::piece<>::rollout(5, axis::x_start(2.0), axis::width_start(1.1), - axis::bin_trans::multiply(1.2)); - BOOST_TEST_EQ(p.size(), 5); - BOOST_TEST(p.bt().is_multiply()); - BOOST_TEST_EQ(p.bt().value(), 1.2); + const auto p_uniform = axis::piece_uniform::solve_b0(8, -1.5, +2.5); + auto p = axis::piecewise>(p_uniform); - BOOST_TEST_EQ(p.b0(), 1.1); - BOOST_TEST_EQ(p.x0(), 2.0); - BOOST_TEST_EQ(p.xN(), p.x0() + p.b0() * (1 + 1.2 + std::pow(1.2, 2) + - std::pow(1.2, 3) + std::pow(1.2, 4))); + // Add 4 bins to the R expanding by 1.1 each time + p.extrapolate_R>(4, 1.1); - BOOST_TEST_EQ(p.size_ic(), 5); - BOOST_TEST_EQ(p.size(), 5); - BOOST_TEST_EQ(p.accumulated_left_shift(), 0); - BOOST_TEST_EQ(p.accumulated_right_shift(), 0); + BOOST_TEST_EQ(p.size(), 12); + BOOST_TEST_EQ(p.x0(), -1.5); - const auto b0_check = axis::piece<>::calc_b0(p.x0(), p.xN(), p.size_ic(), p.bt()); - BOOST_TEST_IS_CLOSE(b0_check.value(), p.b0(), 1e-12); - - BOOST_TEST_EQ(p.forward(2.0), 0.0); - BOOST_TEST_EQ(p.forward(2.0 + 1.1), 1.0); - BOOST_TEST_EQ(p.forward(2.0 + 1.1 * (1 + 1.2)), 2.0); - BOOST_TEST_EQ(p.forward(2.0 + 1.1 * (1 + 1.2 + std::pow(1.2, 2))), 3.0); - - for (int i = -100; i < 100; ++i) { - const auto d = i * 0.021; - BOOST_TEST_IS_CLOSE(p.reverse(p.forward(d)), d, 1e-12); - } + // The bin spacing of the first piece is 0.5, so the first extrapolated piece is 0.5 + // * 1.1 + BOOST_TEST_EQ(p.xN(), 2.5 + 0.5 * (1.1 + std::pow(1.1, 2) + std::pow(1.1, 3) + + std::pow(1.1, 4))); } - // Construct rollout left + // Piecewise Disconnected { - const auto p = axis::piece<>::rollout(-5, axis::x_start(2.0), axis::width_start(1.1), - axis::bin_trans::multiply(1.2)); - BOOST_TEST_EQ(p.size(), 5); - BOOST_TEST(p.bt().is_multiply()); - BOOST_TEST_EQ(p.bt().value(), 1 / 1.2); - - // const auto bt_rev = axis::bin_trans::multiply(1 / 1.2); - BOOST_TEST_IS_CLOSE(p.b0(), 1.1 * std::pow(1.2, 4), 1e-12); - BOOST_TEST_IS_CLOSE(p.xN(), 2.0, 1e-12); - - BOOST_TEST_IS_CLOSE( - p.x0(), - p.xN() - 1.1 * (1 + 1.2 + std::pow(1.2, 2) + std::pow(1.2, 3) + std::pow(1.2, 4)), - 1e-12); - - // const auto bt = p.bt(); - // BOOST_TEST(bt.is_multiply()); - // BOOST_TEST_EQ(bt.value(), 1.2); - - // BOOST_TEST_EQ(p.size_ic(), 5); - // BOOST_TEST_EQ(p.size(), 5); - // BOOST_TEST_EQ(p.accumulated_left_shift(), 0); - // BOOST_TEST_EQ(p.accumulated_right_shift(), 0); - - // std::cout << "p.bt().value(): " << p.bt().value() << std::endl; - // std::cout << "b0: " << p.b0() << std::endl; - // std::cout << "x0: " << p.x0() << std::endl; - // std::cout << "xN: " << p.xN() << std::endl; - - // for (int i = 1; i < 6; ++i) { - // std::cout << "bin_width: " << p.calc_bin_width(i) << std::endl; - // } - - // assert(false); + const auto p_arbitrary = axis::piece_variable::create( + std::vector{-1.5, -1.2, 0.4, 1.2, 2.5}); + auto p = axis::piecewise>(p_arbitrary); + + // New piece + const auto p_uniform = axis::piece_uniform::solve_b0(8, 3.5, +5.5); + p.add_R_gap_okay(p_uniform); + + BOOST_TEST_EQ(p.size(), 12); + BOOST_TEST_EQ(p.x0(), -1.5); + BOOST_TEST_EQ(p.xN(), 5.5); + + // Check indices + BOOST_TEST_EQ(p.forward(-1.5), 0); + BOOST_TEST_EQ(p.forward(-1.2), 1); + BOOST_TEST_EQ(p.forward(0.4), 2); + BOOST_TEST_EQ(p.forward(1.2), 3); + BOOST_TEST_EQ(p.forward(2.49999999999), 3); + + BOOST_TEST(std::isnan(p.forward(2.5))); + BOOST_TEST(std::isnan(p.forward(3.0))); + BOOST_TEST(std::isnan(p.forward(3.49999999999))); + + BOOST_TEST_EQ(p.forward(3.5), 4); + BOOST_TEST_EQ(p.forward(4.5), 8); + BOOST_TEST_EQ(p.forward(5.5), 12); } - // Construct piecewise - { - const auto p = axis::piece<>::rollout(-5, axis::x_start(2.0), axis::width_start(1.1), - axis::bin_trans::multiply(1.2)); - auto pw = axis::piecewise(p); - } - - // // bad_ctors - // { - // BOOST_TEST_THROWS(axis::unit_regular<>(1, 0, 0), std::invalid_argument); - // BOOST_TEST_THROWS(axis::unit_regular<>(0, 0, 1), std::invalid_argument); - // } - - // // ctors and assignment - // { - // axis::unit_regular<> a{4, -2, 2}; - // axis::unit_regular<> b; - // BOOST_TEST_NE(a, b); - // b = a; - // BOOST_TEST_EQ(a, b); - // axis::unit_regular<> c = std::move(b); - // BOOST_TEST_EQ(c, a); - // axis::unit_regular<> d; - // BOOST_TEST_NE(c, d); - // d = std::move(c); - // BOOST_TEST_EQ(d, a); - // } - - // // input, output - // { - // axis::unit_regular<> a{4, -2, 2, "foo"}; - // BOOST_TEST_EQ(a.metadata(), "foo"); - // const auto& cref = a; - // BOOST_TEST_EQ(cref.metadata(), "foo"); - // cref.metadata() = "bar"; // this is allowed - // BOOST_TEST_EQ(cref.metadata(), "bar"); - // BOOST_TEST_EQ(a.value(0), -2); - // BOOST_TEST_EQ(a.value(1), -1); - // BOOST_TEST_EQ(a.value(2), 0); - // BOOST_TEST_EQ(a.value(3), 1); - // BOOST_TEST_EQ(a.value(4), 2); - // BOOST_TEST_EQ(a.bin(-1).lower(), -std::numeric_limits::infinity()); - // BOOST_TEST_EQ(a.bin(-1).upper(), -2); - // BOOST_TEST_EQ(a.bin(a.size()).lower(), 2); - // BOOST_TEST_EQ(a.bin(a.size()).upper(), std::numeric_limits::infinity()); - // BOOST_TEST_EQ(a.index(-10.), -1); - // BOOST_TEST_EQ(a.index(-2.1), -1); - // BOOST_TEST_EQ(a.index(-2.0), 0); - // BOOST_TEST_EQ(a.index(-1.1), 0); - // BOOST_TEST_EQ(a.index(0.0), 2); - // BOOST_TEST_EQ(a.index(0.9), 2); - // BOOST_TEST_EQ(a.index(1.0), 3); - // BOOST_TEST_EQ(a.index(10.), 4); - // BOOST_TEST_EQ(a.index(-std::numeric_limits::infinity()), -1); - // BOOST_TEST_EQ(a.index(std::numeric_limits::infinity()), 4); - // BOOST_TEST_EQ(a.index(std::numeric_limits::quiet_NaN()), 4); - - // BOOST_TEST_EQ(str(a), - // "regular(4, -2, 2, metadata=\"bar\", options=underflow | overflow)"); - // } - - // // with inverted range - // { - // axis::unit_regular<> a{2, 1, -2}; - // BOOST_TEST_EQ(a.bin(-1).lower(), std::numeric_limits::infinity()); - // BOOST_TEST_EQ(a.bin(0).lower(), 1); - // BOOST_TEST_EQ(a.bin(1).lower(), -0.5); - // BOOST_TEST_EQ(a.bin(2).lower(), -2); - // BOOST_TEST_EQ(a.bin(2).upper(), -std::numeric_limits::infinity()); - // BOOST_TEST_EQ(a.index(2), -1); - // BOOST_TEST_EQ(a.index(1.001), -1); - // BOOST_TEST_EQ(a.index(1), 0); - // BOOST_TEST_EQ(a.index(0), 0); - // BOOST_TEST_EQ(a.index(-0.499), 0); - // BOOST_TEST_EQ(a.index(-0.5), 1); - // BOOST_TEST_EQ(a.index(-1), 1); - // BOOST_TEST_EQ(a.index(-2), 2); - // BOOST_TEST_EQ(a.index(-20), 2); - // } - - // // with log transform - // { - // auto a = axis::unit_regular{2, 1e0, 1e2}; - // BOOST_TEST_EQ(a.bin(-1).lower(), 0.0); - // BOOST_TEST_IS_CLOSE(a.bin(0).lower(), 1.0, 1e-9); - // BOOST_TEST_IS_CLOSE(a.bin(1).lower(), 10.0, 1e-9); - // BOOST_TEST_IS_CLOSE(a.bin(2).lower(), 100.0, 1e-9); - // BOOST_TEST_EQ(a.bin(2).upper(), std::numeric_limits::infinity()); - - // BOOST_TEST_EQ(a.index(-1), 2); // produces NaN in conversion - // BOOST_TEST_EQ(a.index(0), -1); - // BOOST_TEST_EQ(a.index(1), 0); - // BOOST_TEST_EQ(a.index(9), 0); - // BOOST_TEST_EQ(a.index(10), 1); - // BOOST_TEST_EQ(a.index(90), 1); - // BOOST_TEST_EQ(a.index(100), 2); - // BOOST_TEST_EQ(a.index(std::numeric_limits::infinity()), 2); - - // BOOST_TEST_THROWS((axis::unit_regular{2, -1, 0}), - // std::invalid_argument); - - // BOOST_TEST_CSTR_EQ( - // str(a).c_str(), - // "regular(transform::log{}, 2, 1, 100, options=underflow | overflow)"); - // } - - // // with sqrt transform - // { - // axis::unit_regular a(2, 0, 4); - // // this is weird, but -inf * -inf = inf, thus the lower bound - // BOOST_TEST_EQ(a.bin(-1).lower(), std::numeric_limits::infinity()); - // BOOST_TEST_IS_CLOSE(a.bin(0).lower(), 0.0, 1e-9); - // BOOST_TEST_IS_CLOSE(a.bin(1).lower(), 1.0, 1e-9); - // BOOST_TEST_IS_CLOSE(a.bin(2).lower(), 4.0, 1e-9); - // BOOST_TEST_EQ(a.bin(2).upper(), std::numeric_limits::infinity()); - - // BOOST_TEST_EQ(a.index(-1), 2); // produces NaN in conversion - // BOOST_TEST_EQ(a.index(0), 0); - // BOOST_TEST_EQ(a.index(0.99), 0); - // BOOST_TEST_EQ(a.index(1), 1); - // BOOST_TEST_EQ(a.index(3.99), 1); - // BOOST_TEST_EQ(a.index(4), 2); - // BOOST_TEST_EQ(a.index(100), 2); - // BOOST_TEST_EQ(a.index(std::numeric_limits::infinity()), 2); - - // BOOST_TEST_EQ(str(a), - // "regular(transform::sqrt{}, 2, 0, 4, options=underflow | overflow)"); - // } - - // // with pow transform - // { - // axis::unit_regular a(tr::pow{0.5}, 2, 0, 4); - // // this is weird, but -inf * -inf = inf, thus the lower bound - // BOOST_TEST_EQ(a.bin(-1).lower(), std::numeric_limits::infinity()); - // BOOST_TEST_IS_CLOSE(a.bin(0).lower(), 0.0, 1e-9); - // BOOST_TEST_IS_CLOSE(a.bin(1).lower(), 1.0, 1e-9); - // BOOST_TEST_IS_CLOSE(a.bin(2).lower(), 4.0, 1e-9); - // BOOST_TEST_EQ(a.bin(2).upper(), std::numeric_limits::infinity()); - - // BOOST_TEST_EQ(a.index(-1), 2); // produces NaN in conversion - // BOOST_TEST_EQ(a.index(0), 0); - // BOOST_TEST_EQ(a.index(0.99), 0); - // BOOST_TEST_EQ(a.index(1), 1); - // BOOST_TEST_EQ(a.index(3.99), 1); - // BOOST_TEST_EQ(a.index(4), 2); - // BOOST_TEST_EQ(a.index(100), 2); - // BOOST_TEST_EQ(a.index(std::numeric_limits::infinity()), 2); - - // BOOST_TEST_EQ(str(a), - // "regular(transform::pow{0.5}, 2, 0, 4, options=underflow | - // overflow)"); - // } - - // // with step - // { - // axis::unit_regular<> a(axis::step(0.5), 1, 3); - // BOOST_TEST_EQ(a.size(), 4); - // BOOST_TEST_EQ(a.bin(-1).lower(), -std::numeric_limits::infinity()); - // BOOST_TEST_EQ(a.value(0), 1); - // BOOST_TEST_EQ(a.value(1), 1.5); - // BOOST_TEST_EQ(a.value(2), 2); - // BOOST_TEST_EQ(a.value(3), 2.5); - // BOOST_TEST_EQ(a.value(4), 3); - // BOOST_TEST_EQ(a.bin(4).upper(), std::numeric_limits::infinity()); - - // axis::unit_regular<> b(axis::step(0.5), 1, 3.1); - // BOOST_TEST_EQ(a, b); - // } - - // // with circular option - // { - // axis::circular<> a{4, 0, 1}; - // BOOST_TEST_EQ(a.bin(-1).lower(), a.bin(a.size() - 1).lower() - 1); - // BOOST_TEST_EQ(a.index(-1.0 * 3), 0); - // BOOST_TEST_EQ(a.index(0.0), 0); - // BOOST_TEST_EQ(a.index(0.25), 1); - // BOOST_TEST_EQ(a.index(0.5), 2); - // BOOST_TEST_EQ(a.index(0.75), 3); - // BOOST_TEST_EQ(a.index(1.0), 0); - // BOOST_TEST_EQ(a.index(std::numeric_limits::infinity()), 4); - // BOOST_TEST_EQ(a.index(-std::numeric_limits::infinity()), 4); - // BOOST_TEST_EQ(a.index(std::numeric_limits::quiet_NaN()), 4); - // } - - // // with growth - // { - // using pii_t = std::pair; - // axis::unit_regular a{1, 0, 1}; - // BOOST_TEST_EQ(a.size(), 1); - // BOOST_TEST_EQ(a.update(0), pii_t(0, 0)); - // BOOST_TEST_EQ(a.size(), 1); - // BOOST_TEST_EQ(a.update(1), pii_t(1, -1)); - // BOOST_TEST_EQ(a.size(), 2); - // BOOST_TEST_EQ(a.value(0), 0); - // BOOST_TEST_EQ(a.value(2), 2); - // BOOST_TEST_EQ(a.update(-1), pii_t(0, 1)); - // BOOST_TEST_EQ(a.size(), 3); - // BOOST_TEST_EQ(a.value(0), -1); - // BOOST_TEST_EQ(a.value(3), 2); - // BOOST_TEST_EQ(a.update(-10), pii_t(0, 9)); - // BOOST_TEST_EQ(a.size(), 12); - // BOOST_TEST_EQ(a.value(0), -10); - // BOOST_TEST_EQ(a.value(12), 2); - // BOOST_TEST_EQ(a.update(std::numeric_limits::infinity()), pii_t(a.size(), - // 0)); BOOST_TEST_EQ(a.update(std::numeric_limits::quiet_NaN()), - // pii_t(a.size(), 0)); - // BOOST_TEST_EQ(a.update(-std::numeric_limits::infinity()), pii_t(-1, 0)); - // } - - // // axis with overflow bin represents open interval - // { - // axis::unit_regular a{2, 0, 1}; - // BOOST_TEST_EQ(a.index(0), 0); - // BOOST_TEST_EQ(a.index(0.49), 0); - // BOOST_TEST_EQ(a.index(0.50), 1); - // BOOST_TEST_EQ(a.index(0.99), 1); - // BOOST_TEST_EQ(a.index(1), 2); // overflow bin - // BOOST_TEST_EQ(a.index(1.1), 2); // overflow bin - // } - - // // axis without overflow bin represents a closed interval - // { - // axis::unit_regular a{2, 0, 1}; - // BOOST_TEST_EQ(a.index(0), 0); - // BOOST_TEST_EQ(a.index(0.49), 0); - // BOOST_TEST_EQ(a.index(0.50), 1); - // BOOST_TEST_EQ(a.index(0.99), 1); - // BOOST_TEST_EQ(a.index(1), 1); // last ordinary bin - // BOOST_TEST_EQ(a.index(1.1), 2); // out of range - // } - - // // iterators - // { - // test_axis_iterator(axis::unit_regular<>(5, 0, 1), 0, 5); - // test_axis_iterator(axis::unit_regular(5, 0, 1), 0, - // 5); test_axis_iterator(axis::circular<>(5, 0, 1), 0, 5); - // } - - // // bin_type streamable - // { - // auto test = [](const auto& x, const char* ref) { - // std::ostringstream os; - // os << x; - // BOOST_TEST_EQ(os.str(), std::string(ref)); - // }; - - // auto a = axis::unit_regular<>(2, 0, 1); - // test(a.bin(0), "[0, 0.5)"); - // } - - // // null_type streamable - // { - // auto a = axis::unit_regular(2, 0, 1); - // BOOST_TEST_EQ(str(a), "regular(2, 0, 1, options=underflow | overflow)"); - // } - - // // shrink and rebin - // { - // using A = axis::unit_regular<>; - // auto a = A(5, 0, 5); - // auto b = A(a, 1, 4, 1); - // BOOST_TEST_EQ(b.size(), 3); - // BOOST_TEST_EQ(b.value(0), 1); - // BOOST_TEST_EQ(b.value(3), 4); - // auto c = A(a, 0, 4, 2); - // BOOST_TEST_EQ(c.size(), 2); - // BOOST_TEST_EQ(c.value(0), 0); - // BOOST_TEST_EQ(c.value(2), 4); - // auto e = A(a, 1, 5, 2); - // BOOST_TEST_EQ(e.size(), 2); - // BOOST_TEST_EQ(e.value(0), 1); - // BOOST_TEST_EQ(e.value(2), 5); - // } - - // // shrink and rebin with circular option - // { - // using A = axis::circular<>; - // auto a = A(4, 1, 5); - // BOOST_TEST_THROWS(A(a, 1, 4, 1), std::invalid_argument); - // BOOST_TEST_THROWS(A(a, 0, 3, 1), std::invalid_argument); - // auto b = A(a, 0, 4, 2); - // BOOST_TEST_EQ(b.size(), 2); - // BOOST_TEST_EQ(b.value(0), 1); - // BOOST_TEST_EQ(b.value(2), 5); - // } - return boost::report_errors(); } From f0625b4440d9f74b668a79a55e858b3de4efa378 Mon Sep 17 00:00:00 2001 From: Ryan Date: Tue, 13 Jun 2023 15:25:42 -0400 Subject: [PATCH 04/10] removed overview and select unneeded headers --- include/boost/histogram/axis/int_resolver.hpp | 1 - include/boost/histogram/axis/piece.hpp | 3 +- include/boost/histogram/axis/piecewise.hpp | 89 ------------------- test/axis_int_resolver_test.cpp | 9 +- test/axis_piece_test.cpp | 8 +- test/axis_piecewise_test.cpp | 7 +- 6 files changed, 9 insertions(+), 108 deletions(-) diff --git a/include/boost/histogram/axis/int_resolver.hpp b/include/boost/histogram/axis/int_resolver.hpp index a9b7e130..00d0f4b8 100644 --- a/include/boost/histogram/axis/int_resolver.hpp +++ b/include/boost/histogram/axis/int_resolver.hpp @@ -23,7 +23,6 @@ #include #include #include -#include #include #include #include diff --git a/include/boost/histogram/axis/piece.hpp b/include/boost/histogram/axis/piece.hpp index 92098bb9..61f68aad 100644 --- a/include/boost/histogram/axis/piece.hpp +++ b/include/boost/histogram/axis/piece.hpp @@ -13,7 +13,7 @@ #include #include #include -#include +#include // For one_unit, get_scale_type, etc #include #include #include @@ -33,7 +33,6 @@ namespace boost { namespace histogram { namespace axis { - /** The abstract piece class that other pieces inherit from The class has two pure virtual functions: diff --git a/include/boost/histogram/axis/piecewise.hpp b/include/boost/histogram/axis/piecewise.hpp index a131154a..89e583d9 100644 --- a/include/boost/histogram/axis/piecewise.hpp +++ b/include/boost/histogram/axis/piecewise.hpp @@ -22,7 +22,6 @@ #include #include #include -#include #include #include #include @@ -35,94 +34,6 @@ namespace boost { namespace histogram { namespace axis { -/** - -Solution overview: - -This 1D piecewise axis implementation has two layers. The first layer has to -do with creating immutable transformations from input space X to bin space Y. -The second layer modifies the behavior of these transformation in order to make -them behave like an axis. If axis growth is allowed, this second layer is mutable. -I explain how this works on a high level. - -Layer 1: Transformation -Transformation pieces: - -Suppose we want to divide the input space X into 4 bins of unit spacing starting -at x=3. - -bin 0 1 2 3 - |---|---|---|---| - x=3 x=7 - -The transformation is described by three pieces of information - forward transformation: y = f(x) = x - 3 - inverse transformation: x = f⁻¹(y) = y + 3 - number of bins: N = 4 - -There are five supported transformation types: - 1. unit (bin spacing is 1) - 2. uniform (bin spacing is constant) - 3. multiply (bin spacing is multiplied by a constant) - 4. add (bin spacing is added to by a constant) - 5. arbitrary (bin bounds are user supplied points) - - -Combining transformations: - -Suppose we wanted to combine the following two transformations. - - bin 0 1 2 3 bin 0 1 2 3 4 - |---|---|---|---| |-|--|---|----|-----| - x=3 x=7 x=9 x=15 - Unit spacing Add spacing - -Typically one wouldn't want to combine two transformations with a gap between them, -but this possibility is supported by the implementation. A piecewise transformation -stores each of these unmodified transformations internally. The external behavior -is that of transformation with bins 0 to 8 shown below. - - bin 0 1 2 3 nan 4 5 6 7 8 - |---|---|---|---| |-|--|---|----|-----| - x=3 x=7 x=9 x=15 - - -Layer 2: Axis behavior -Bin shifting: - -Suppose we wanted to use the transformation above in a growable axis. This axis starts -off with 9 bins. Two kinds of bins can be added: underflow and overflow. Underflow bins -are added below x=3. Overflow bins are added above x=15. Adding bins between x=7 and -x=9 is not supported. The mutable transformation `x_bin_shift` accumulates: - - added underflow bins, and - - added overflow bins. -Using this information, x_bin_shift augments the behavior of the contained piecewise -transformation as follows. The size (i.e., number of bins) of the transformation is -increased by the number of underflow and overflow bins. The bin values (i.e., Y axis) -is shifted by the number of underflow bins in the forward and inverse calculations. - -Linear axis: - -Wrapping the, possibly shifted, transformation in an `int_resolver_linear` makes the -transformation behave like an axis. This linear int resolver implemenets an index method. - -Circular axis: - -Suppose we wanted to make the above described piecewise transformation circular with -a relevant period x=3 to x=23. This transformaion is copied below for convenience. - - bin 0 1 2 3 nan 4 5 6 7 8 - |---|---|---|---| |-|--|---|----|-----| - x=3 x=7 x=9 x=15 - -To do this, we wrap the transformation in an `int_resolver_circular` with the bounds -x=3 and x=23. This circular transform augments the behavior of the contained -transformation. The input (i.e., the X axis) is wrapped to the domain x=3 to x=23 -before being given to the contained transformation. This circular int resolver -implemenets an index method allowing it to act as an axis. - -*/ - /** Piecewise -- putting pieces together */ diff --git a/test/axis_int_resolver_test.cpp b/test/axis_int_resolver_test.cpp index d40b19c4..adfc18c6 100644 --- a/test/axis_int_resolver_test.cpp +++ b/test/axis_int_resolver_test.cpp @@ -19,7 +19,6 @@ #include "throw_exception.hpp" #include -#include #include using namespace boost::histogram; @@ -29,7 +28,7 @@ int main() { { using type_multiply = axis::piece_multiply; - // Has x0 = 4, x1 = 5, x2 = 7, x3 = 11, x4 = 17 + // Has x0 = 3, x1 = 4, x2 = 6, x3 = 10, x4 = 18 const int k_N = 3; const double k_x0 = 3.0; const double k_b0 = 1.0; @@ -43,7 +42,6 @@ int main() { BOOST_TEST_EQ(a.index(6), 2); BOOST_TEST_EQ(a.index(10), 3); BOOST_TEST_EQ(a.index(18), 4); - // BOOST_TEST_EQ(a.index(5), k_x0 + k_b0); } // int_resolver_linear @@ -58,10 +56,7 @@ int main() { BOOST_TEST_EQ(a.value(2), 0); BOOST_TEST_EQ(a.value(3), 1); BOOST_TEST_EQ(a.value(4), 2); - // BOOST_TEST_EQ(a.bin(-1).lower(), -std::numeric_limits::infinity()); - // BOOST_TEST_EQ(a.bin(-1).upper(), -2); - // BOOST_TEST_EQ(a.bin(a.size()).lower(), 2); - // BOOST_TEST_EQ(a.bin(a.size()).upper(), std::numeric_limits::infinity()); + BOOST_TEST_EQ(a.index(-10.), -1); BOOST_TEST_EQ(a.index(-2.1), -1); BOOST_TEST_EQ(a.index(-2.0), 0); diff --git a/test/axis_piece_test.cpp b/test/axis_piece_test.cpp index 0fecb153..1fe649d4 100644 --- a/test/axis_piece_test.cpp +++ b/test/axis_piece_test.cpp @@ -17,7 +17,6 @@ #include "throw_exception.hpp" #include -#include #include using namespace boost::histogram; @@ -52,7 +51,6 @@ class tester_piece_b0 { // Check that the bin width is correct BOOST_TEST_IS_CLOSE(p.calc_bin_width(i), bin_width, 1.0e-12); bin_width = lambda(bin_width); - // tester_next_width(p, &bin_width, args...); } // Check xN @@ -74,6 +72,7 @@ class tester_piece_b0 { double b0_; }; +// Tests that two pieces are equal template void test_equal(const PieceType& p1, const PieceType& p2) { BOOST_TEST_EQ(p1.size(), p2.size()); @@ -88,9 +87,10 @@ void test_equal(const PieceType& p1, const PieceType& p2) { template void test_solve_b0(const PieceType& p_R, Args... args) { - // Construct force from create + // Reconstruct piece with xN instead of b0 const auto p2 = PieceType::solve_b0(p_R.size(), p_R.x0(), p_R.xN(), args...); + // See if pieces are equal test_equal(p_R, p2); } @@ -120,7 +120,7 @@ int main() { BOOST_TEST_EQ(p.inverse(4), 6); } - // Test arbitrary + // Test variable { using p_type = axis::piece_variable; const auto p = p_type::create(std::vector{1.1, 2.6, 3.8, 4.2, 5.9}); diff --git a/test/axis_piecewise_test.cpp b/test/axis_piecewise_test.cpp index bf3d73f6..0af0bb65 100644 --- a/test/axis_piecewise_test.cpp +++ b/test/axis_piecewise_test.cpp @@ -18,11 +18,8 @@ #include "throw_exception.hpp" #include -#include #include -#include - using namespace boost::histogram; int main() { @@ -37,8 +34,8 @@ int main() { BOOST_TEST_EQ(p.size(), 12); BOOST_TEST_EQ(p.x0(), -1.5); - // The bin spacing of the first piece is 0.5, so the first extrapolated piece is 0.5 - // * 1.1 + // The bin spacing of the first piece is 0.5, so the first extrapolated piece is + // 0.5 * 1.1 BOOST_TEST_EQ(p.xN(), 2.5 + 0.5 * (1.1 + std::pow(1.1, 2) + std::pow(1.1, 3) + std::pow(1.1, 4))); } From eea42801458146ceccf16406e1c349297334444a Mon Sep 17 00:00:00 2001 From: Ryan Date: Tue, 13 Jun 2023 15:48:34 -0400 Subject: [PATCH 05/10] minor documentation improvements --- include/boost/histogram/axis/int_resolver.hpp | 13 ++++++------ include/boost/histogram/axis/piece.hpp | 21 ++++++++----------- include/boost/histogram/fwd.hpp | 6 ------ 3 files changed, 15 insertions(+), 25 deletions(-) diff --git a/include/boost/histogram/axis/int_resolver.hpp b/include/boost/histogram/axis/int_resolver.hpp index 00d0f4b8..50258c5e 100644 --- a/include/boost/histogram/axis/int_resolver.hpp +++ b/include/boost/histogram/axis/int_resolver.hpp @@ -37,7 +37,8 @@ namespace axis { /** Bin shift integrator - Accumulates bin shifts to support growable axes. + Accumulates bin shifts to support growable axes. The shift_axis method is called + by the update method of an int_resolver class. */ template class bin_shift_integrator { @@ -48,7 +49,7 @@ class bin_shift_integrator { /// Constructor for bin_shift_integrator explicit bin_shift_integrator(const Payload& payload) : payload_(payload) {} - /// Shifts the axis + /// Shifts the axis -- called by update method of int_resolver void shift_axis(index_type n) { if (n < 0) { bins_under_ += std::abs(n); @@ -94,8 +95,6 @@ class int_resolver_linear { /// The mapping from input space X to integer bin space Y index_type index(value_type x) const noexcept { - // Runs in hot loop, please measure impact of changes - const value_type y = payload_.forward(x); if (y < size()) { @@ -185,16 +184,16 @@ class int_resolver_circular { } index_type index(value_type x) const noexcept { - // Take the mod of the input + // Take the mod of the input. Taking the mod of the output likely doesn't do what + // one wants for non-linear transformations. x -= x_low_plus_delta_x_ * std::floor((x - x_low_) / x_low_plus_delta_x_); + // x is in undefined part of periodic domain --> nan if (x < x_low_ || x_high_ < x) { return payload_.size(); } else { const value_type y = payload_.forward(x); - if (std::isfinite(y)) { return static_cast(y); } - return payload_.size(); } } diff --git a/include/boost/histogram/axis/piece.hpp b/include/boost/histogram/axis/piece.hpp index 61f68aad..90fcbb57 100644 --- a/include/boost/histogram/axis/piece.hpp +++ b/include/boost/histogram/axis/piece.hpp @@ -167,8 +167,10 @@ class piece_variable : public piece_base> { std::vector vec_; }; -/** Abstract variable width piece class +/** Abstract class for b0 pieces + Pieces derived from this class have forward and reverse methods with closed form + solutions parameterized in terms of an initial bin spacing b0. */ template class piece_b0 : public piece_base { @@ -176,7 +178,7 @@ class piece_b0 : public piece_base { using internal_value_type = detail::get_scale_type; public: - /** Finds the piece whose size is known, but bin spacing is not. + /** Finds the piece whose size is known, but initial bin spacing b0 is not. @param n number of bins @param x0 starting point @@ -289,14 +291,6 @@ class piece_uniform : public piece_b0> { x1 - x0 = r * b0 / (r - 1) - 1 * b0 / (r - 1) = (r - 1) * b0 / (r - 1) = b0 - - Find a formula for N - xN - x0 = b0 / (r - 1) * rᴺ - b0 / (r - 1) * 1 - = b0 / (r - 1) * (rᴺ - 1) - - N in terms of xN - N = log2(((xN - x0) / b0) * (r - 1) + 1) / log2(r) - */ template class piece_multiply : public piece_b0> { @@ -307,8 +301,11 @@ class piece_multiply : public piece_b0> { public: /** The mapping from bin space Y to input space X. - Uses the formula: - x = x0 + b0 * (rᴺ - 1) / (r - 1) + Eliminate ψ: + xN - x0 = b0 / (r - 1) * rᴺ - b0 / (r - 1) * 1 + = b0 / (r - 1) * (rᴺ - 1) + Solving for xN gives the formula: + xN = x0 + b0 * (rᴺ - 1) / (r - 1) */ internal_value_type inverse(internal_value_type N) const noexcept override { return this->x0() + this->b0() * (std::pow(r_, N) - 1) / (r_ - 1); diff --git a/include/boost/histogram/fwd.hpp b/include/boost/histogram/fwd.hpp index 3340840f..3c48f784 100644 --- a/include/boost/histogram/fwd.hpp +++ b/include/boost/histogram/fwd.hpp @@ -60,12 +60,6 @@ template class regular; -template -class piece; - -template -class unit_regular; - template class integer; From c59ee4dc5e43d1ea2479a9a0efaee1de7aad2890 Mon Sep 17 00:00:00 2001 From: Ryan Date: Tue, 13 Jun 2023 15:54:23 -0400 Subject: [PATCH 06/10] added wraparound to int_resolver test --- test/axis_int_resolver_test.cpp | 13 ++++++++----- 1 file changed, 8 insertions(+), 5 deletions(-) diff --git a/test/axis_int_resolver_test.cpp b/test/axis_int_resolver_test.cpp index adfc18c6..eed137a1 100644 --- a/test/axis_int_resolver_test.cpp +++ b/test/axis_int_resolver_test.cpp @@ -37,11 +37,14 @@ int main() { const auto p = type_multiply::create(k_N, k_x0, k_b0, k_r); const auto a = axis::int_resolver_circular(p, 2, 22); - BOOST_TEST_EQ(a.index(3), 0); - BOOST_TEST_EQ(a.index(4), 1); - BOOST_TEST_EQ(a.index(6), 2); - BOOST_TEST_EQ(a.index(10), 3); - BOOST_TEST_EQ(a.index(18), 4); + // Has a period of 20 based on the given bounds 2 and 22 + for (int i = -200; i < 200; i += 20) { + BOOST_TEST_EQ(a.index(i + 3), 0); + BOOST_TEST_EQ(a.index(i + 4), 1); + BOOST_TEST_EQ(a.index(i + 6), 2); + BOOST_TEST_EQ(a.index(i + 10), 3); + BOOST_TEST_EQ(a.index(i + 18), 4); + } } // int_resolver_linear From 402ab3dfafe8f6ed48f0d0263a8697cad1dd8c65 Mon Sep 17 00:00:00 2001 From: Ryan Date: Tue, 13 Jun 2023 17:12:27 -0400 Subject: [PATCH 07/10] See if auto fixes compile issue on AppVeyor --- include/boost/histogram/axis/piecewise.hpp | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/include/boost/histogram/axis/piecewise.hpp b/include/boost/histogram/axis/piecewise.hpp index 89e583d9..051ae924 100644 --- a/include/boost/histogram/axis/piecewise.hpp +++ b/include/boost/histogram/axis/piecewise.hpp @@ -51,7 +51,7 @@ class piecewise { /// The mapping from input space X to bin space Y template - T forward(T x) const noexcept { + auto forward(T x) const noexcept { // First piece const auto& piece_0 = v_pieces_.front(); From aac64f142eb1636e2e026b940032f5174abe8656 Mon Sep 17 00:00:00 2001 From: Ryan Date: Wed, 14 Jun 2023 13:30:31 -0400 Subject: [PATCH 08/10] Added test based on example case in #336 --- test/axis_int_resolver_test.cpp | 44 +++++++++++++++++++++++++++++++++ 1 file changed, 44 insertions(+) diff --git a/test/axis_int_resolver_test.cpp b/test/axis_int_resolver_test.cpp index eed137a1..9cc57690 100644 --- a/test/axis_int_resolver_test.cpp +++ b/test/axis_int_resolver_test.cpp @@ -105,5 +105,49 @@ int main() { BOOST_TEST_EQ(a.update(-std::numeric_limits::infinity()), pii_t(-1, 0)); } + // Test based on the example https://godbolt.org/z/oaPo6n17h + // from @vakokako in #336 + { + auto fn_test_precision = [](int N, double x0, double xN, auto fn_axis) { + const auto a = fn_axis(N, x0, xN); + + // // Calculate bin spacing b0 + const double b0 = (xN - x0) / N; + + // Check to see if the index and value calculations are exact + for (int y = 0; y < a.size(); ++y) { + const double x = x0 + y * b0; + BOOST_TEST_EQ(y, a.index(x)); + BOOST_TEST_EQ((double)(x), a.value(y)); + } + }; + + auto fn_test_piece = [](int N, double x0, double xN) { + const auto p_uniform = axis::piece_uniform::solve_b0(N, x0, xN); + return axis::int_resolver_linear>(p_uniform); + }; + +#if 0 + auto fn_test_regular = [](int N, double x0, double xN){ + return boost::histogram::axis::regular(N, x0, xN); + }; + fn_test_precision(27000, 0, 27000, fn_test_regular); // Fails +#endif + fn_test_precision(27000, 0, 27000, fn_test_piece); // Passes + + // Bin spacings and starting points that take few floating point bits to represent + const std::vector v_spacing = {0.125, 0.25, 0.375, 0.5, 0.75, 1, 1.5, 3, 7.5}; + const std::vector v_x0 = {-1000.25, -2.5, -0.5, 0, 0.5, 2.5, 1000.25}; + const std::vector v_n = {1, 16, 27000, 350017, 1234567}; // Integer spacings + + for (const double spacing : v_spacing) { + for (const double x0 : v_x0) { + for (const int n : v_n) { + fn_test_precision(n, x0, x0 + n * spacing, fn_test_piece); + } + } + } + } + return boost::report_errors(); } From 8b2ff882c7c199441133dcfc5984e5f1d2c733f8 Mon Sep 17 00:00:00 2001 From: Ryan Date: Wed, 14 Jun 2023 13:40:02 -0400 Subject: [PATCH 09/10] Fix incorrect comment --- test/axis_int_resolver_test.cpp | 9 ++++----- 1 file changed, 4 insertions(+), 5 deletions(-) diff --git a/test/axis_int_resolver_test.cpp b/test/axis_int_resolver_test.cpp index 9cc57690..a2351e71 100644 --- a/test/axis_int_resolver_test.cpp +++ b/test/axis_int_resolver_test.cpp @@ -136,13 +136,12 @@ int main() { fn_test_precision(27000, 0, 27000, fn_test_piece); // Passes // Bin spacings and starting points that take few floating point bits to represent - const std::vector v_spacing = {0.125, 0.25, 0.375, 0.5, 0.75, 1, 1.5, 3, 7.5}; + const std::vector v_spacing = {0.125, 0.25, 0.375, 0.5, 0.75, 1, 1.625, 7.25}; const std::vector v_x0 = {-1000.25, -2.5, -0.5, 0, 0.5, 2.5, 1000.25}; - const std::vector v_n = {1, 16, 27000, 350017, 1234567}; // Integer spacings - for (const double spacing : v_spacing) { - for (const double x0 : v_x0) { - for (const int n : v_n) { + for (const int n : {1, 16, 27000, 350017, 1234567}) { + for (const double spacing : v_spacing) { + for (const double x0 : v_x0) { fn_test_precision(n, x0, x0 + n * spacing, fn_test_piece); } } From 28838d2cf50982ed7546da9e86b77d9d810275b2 Mon Sep 17 00:00:00 2001 From: Ryan Date: Wed, 14 Jun 2023 14:11:25 -0400 Subject: [PATCH 10/10] added size sanity check --- test/axis_int_resolver_test.cpp | 3 ++- 1 file changed, 2 insertions(+), 1 deletion(-) diff --git a/test/axis_int_resolver_test.cpp b/test/axis_int_resolver_test.cpp index a2351e71..db1c30a9 100644 --- a/test/axis_int_resolver_test.cpp +++ b/test/axis_int_resolver_test.cpp @@ -110,6 +110,7 @@ int main() { { auto fn_test_precision = [](int N, double x0, double xN, auto fn_axis) { const auto a = fn_axis(N, x0, xN); + BOOST_TEST(a.size() == N); // // Calculate bin spacing b0 const double b0 = (xN - x0) / N; @@ -136,7 +137,7 @@ int main() { fn_test_precision(27000, 0, 27000, fn_test_piece); // Passes // Bin spacings and starting points that take few floating point bits to represent - const std::vector v_spacing = {0.125, 0.25, 0.375, 0.5, 0.75, 1, 1.625, 7.25}; + const std::vector v_spacing = {0.125, 0.375, 0.5, 0.75, 1, 1.625, 3, 7.25}; const std::vector v_x0 = {-1000.25, -2.5, -0.5, 0, 0.5, 2.5, 1000.25}; for (const int n : {1, 16, 27000, 350017, 1234567}) {