diff --git a/include/boost/histogram/axis.hpp b/include/boost/histogram/axis.hpp index ae42582b..bc95f80c 100644 --- a/include/boost/histogram/axis.hpp +++ b/include/boost/histogram/axis.hpp @@ -19,7 +19,10 @@ #include #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..50258c5e --- /dev/null +++ b/include/boost/histogram/axis/int_resolver.hpp @@ -0,0 +1,212 @@ +// 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 + +namespace boost { +namespace histogram { +namespace axis { + +/** Bin shift integrator + + 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 { + 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 -- called by update method of int_resolver + 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 { + 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. 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(); + } + } + +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..90fcbb57 --- /dev/null +++ b/include/boost/histogram/axis/piece.hpp @@ -0,0 +1,508 @@ +// 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 // For one_unit, get_scale_type, etc +#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 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 { + using value_type = Value; + using internal_value_type = detail::get_scale_type; + +public: + /** Finds the piece whose size is known, but initial bin spacing b0 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 +*/ +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. + + 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); + } + + /** 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 new file mode 100644 index 00000000..051ae924 --- /dev/null +++ b/include/boost/histogram/axis/piecewise.hpp @@ -0,0 +1,164 @@ +// 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 +#include +#include +#include +#include + +namespace boost { +namespace histogram { +namespace axis { + +/** Piecewise -- putting pieces together + +*/ +template +class piecewise { + using value_type = Value; + using internal_value_type = detail::get_scale_type; + +public: + /// 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 + auto forward(T x) const noexcept { + + // First piece + const auto& piece_0 = v_pieces_.front(); + if (x < piece_0.xN()) return piece_0.forward(x); + + // 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]; + + if (x < p.x0()) { + return std::numeric_limits::quiet_NaN(); // Between discontinuous pieces + } else if (x < p.xN()) { + return offset + p.forward(x); + } + + offset += p.size(); + } + + const auto& piece_N = get_piece_R(); + if (piece_N.x0() <= x) { + return offset + get_piece_R().forward(x); // Last piece + } else { + return std::numeric_limits::quiet_NaN(); // Before last piece + } + } + + /** Extrapolates the bin spacing of the right side of the piecewise axis. + + @param N number of bins to extrapolate + @param args arguments to pass to the piece constructor + + 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. + + this->extrapolate_R(3, 2.0); + + <- - - -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); + } + + /** Adds a piece to the right side of the piecewise axis. + + @param p piece to add + + 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); + } + + /** Adds a piece to the right side of the piecewise axis. + + @param p piece to add + + 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")); + + 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; + } + + /// The location of the right side of the right piece. + double xN() const noexcept { return get_piece_R().xN(); } + + /// The location of the left side of the left piece. + double x0() const noexcept { return v_pieces_.front().x0(); } + +private: + // Gets the right-most piece. + const PieceType& get_piece_R() const noexcept { + assert(!v_pieces_.empty()); + return v_pieces_.back(); + } + + std::vector v_pieces_{}; +}; + +} // namespace axis +} // namespace histogram +} // namespace boost + +#endif diff --git a/test/CMakeLists.txt b/test/CMakeLists.txt index 337da698..f208a82a 100644 --- a/test/CMakeLists.txt +++ b/test/CMakeLists.txt @@ -54,6 +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 d38b5398..7523d39e 100644 --- a/test/Jamfile +++ b/test/Jamfile @@ -57,6 +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..db1c30a9 --- /dev/null +++ b/test/axis_int_resolver_test.cpp @@ -0,0 +1,153 @@ +// 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 + +using namespace boost::histogram; + +int main() { + // int_resolver_circular + { + using type_multiply = axis::piece_multiply; + + // 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; + 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); + + // 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 + { + 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.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)); + } + + // 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); + BOOST_TEST(a.size() == N); + + // // 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.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}) { + for (const double spacing : v_spacing) { + for (const double x0 : v_x0) { + fn_test_precision(n, x0, x0 + n * spacing, fn_test_piece); + } + } + } + } + + return boost::report_errors(); +} diff --git a/test/axis_piece_test.cpp b/test/axis_piece_test.cpp new file mode 100644 index 00000000..1fe649d4 --- /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 + +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); + } + + // 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_; +}; + +// Tests that two pieces are equal +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) { + // 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); +} + +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 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}); + + 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 new file mode 100644 index 00000000..0af0bb65 --- /dev/null +++ b/test/axis_piecewise_test.cpp @@ -0,0 +1,74 @@ +// 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 "axis.hpp" +#include "is_close.hpp" +#include "std_ostream.hpp" +#include "str.hpp" +#include "throw_exception.hpp" + +#include +#include + +using namespace boost::histogram; + +int main() { + // Piecewise Connected + { + const auto p_uniform = axis::piece_uniform::solve_b0(8, -1.5, +2.5); + auto p = axis::piecewise>(p_uniform); + + // Add 4 bins to the R expanding by 1.1 each time + p.extrapolate_R>(4, 1.1); + + 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 + 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))); + } + + // Piecewise Disconnected + { + 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); + } + + return boost::report_errors(); +}