From b2ac83cd1e805c890cbe51b22154c95d14a757c9 Mon Sep 17 00:00:00 2001 From: Ryan Date: Thu, 15 Jun 2023 18:59:13 -0400 Subject: [PATCH 1/5] refactored regular --- include/boost/histogram/axis/regular.hpp | 120 ++++++++++++++--------- 1 file changed, 71 insertions(+), 49 deletions(-) diff --git a/include/boost/histogram/axis/regular.hpp b/include/boost/histogram/axis/regular.hpp index a5c7c124..bfcefa40 100644 --- a/include/boost/histogram/axis/regular.hpp +++ b/include/boost/histogram/axis/regular.hpp @@ -83,18 +83,19 @@ struct id { void serialize(Archive&, unsigned /* version */) {} }; -/// Log transform for equidistant bins in log-space. +/// Log transform for equidistant bins in log-space. Note, the base does not matter +/// because of the change of base formula which says: log_b(x) = log_a(x) / log_a(b). struct log { /// Returns log(x) of external value x. template static T forward(T x) { - return std::log(x); + return std::log2(x); } /// Returns exp(x) for internal value x. template static T inverse(T x) { - return std::exp(x); + return std::pow(2, x); } template @@ -221,8 +222,13 @@ class regular : public iterator_mixin(n)) - , min_(this->forward(detail::get_scale(start))) - , delta_(this->forward(detail::get_scale(stop)) - min_) { + , min_(this->forward(detail::get_scale(start))) { + { + const auto delta = this->forward(detail::get_scale(stop)) - min_; + b0_ = static_cast(delta / size()); + b0_inv_ = static_cast(size() / delta); + } + // static_asserts were moved here from class scope to satisfy deduction in gcc>=11 static_assert(std::is_nothrow_move_constructible::value, "transform must be no-throw move constructible"); @@ -233,11 +239,10 @@ class regular : public iterator_mixin 0 required")); - if (!std::isfinite(min_) || !std::isfinite(delta_)) + if (b0_ == 0) BOOST_THROW_EXCEPTION(std::invalid_argument("range of axis is zero")); + if (!std::isfinite(min_) || !std::isfinite(b0_) || !std::isfinite(b0_inv_)) 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")); } /** Construct n bins over real range [start, stop). @@ -306,52 +311,32 @@ class regular : public iterator_mixinforward(x / unit_type{}) - min_) / delta_; - if (options_type::test(option::circular)) { - if (std::isfinite(z)) { - z -= std::floor(z); - return static_cast(z * size()); - } - } else { - if (z < 1) { - if (z >= 0) - return static_cast(z * size()); - else - return -1; - } - // 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 + return index_impl(options_type::test(axis::option::circular), + static_cast(x / unit_type{})); } /// 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 z = (this->forward(x / unit_type{}) - min_) / delta_; - if (z < 1) { // don't use i here! - if (z >= 0) { - const auto i = static_cast(z * size()); + auto y = (this->forward(x / unit_type{}) - min_) * b0_inv_; + if (y < size()) { + if (0 <= y) { + const auto i = static_cast(y); return {i, 0}; } - if (z != -std::numeric_limits::infinity()) { - const auto stop = min_ + delta_; - const auto i = static_cast(std::floor(z * size())); - min_ += i * (delta_ / size()); - delta_ = stop - min_; + if (y != -std::numeric_limits::infinity()) { + const auto i = static_cast(std::floor(y)); + min_ += i * b0_; size_ -= i; return {0, -i}; } // z is -infinity return {-1, 0}; } - // z either beyond range, infinite, or NaN - if (z < std::numeric_limits::infinity()) { - const auto i = static_cast(z * size()); + // y either beyond range, infinite, or NaN + if (y < std::numeric_limits::infinity()) { + const auto i = static_cast(y); const auto n = i - size() + 1; - delta_ /= size(); - delta_ *= size() + n; size_ += n; return {i, -n}; } @@ -361,13 +346,13 @@ class regular : public iterator_mixin::infinity() * delta_; - else if (options_type::test(option::circular) || z <= 1.0) - z = (1.0 - z) * min_ + z * (min_ + delta_); + real_index_type z; + if (!options_type::test(option::circular) && i < 0.0) + z = -std::numeric_limits::infinity() * b0_; + else if (options_type::test(option::circular) || i <= size()) + z = min_ + i * b0_; else { - z = std::numeric_limits::infinity() * delta_; + z = std::numeric_limits::infinity() * b0_; } return static_cast(this->inverse(z) * unit_type()); } @@ -386,7 +371,7 @@ class regular : public iterator_mixin bool operator==(const regular& o) const noexcept { return detail::relaxed_equal{}(transform(), o.transform()) && size() == o.size() && - min_ == o.min_ && delta_ == o.delta_ && + min_ == o.min_ && b0_ == o.b0_ && b0_inv_ == o.b0_inv_ && detail::relaxed_equal{}(this->metadata(), o.metadata()); } template @@ -400,12 +385,49 @@ class regular : public iterator_mixinmetadata()); ar& make_nvp("min", min_); - ar& make_nvp("delta", delta_); + ar& make_nvp("b0", b0_); + ar& make_nvp("b0_inv", b0_inv_); } private: + // axis not circular + index_type index_impl(std::false_type, double x) const noexcept { + auto y = (this->forward(x) - min_) * b0_inv_; + if (y < size()) { + if (0 <= y) + return static_cast(y); + else + return -1; + } + // upper edge of last bin is inclusive if overflow bin is not present + if (!options_type::test(option::overflow) && y == size()) return size() - 1; + return size(); // also returned if x is NaN + } + + // value_type is integer, axis circular + index_type index_impl(std::true_type, double x) const noexcept { + auto lambda_circ = [&](auto x) -> index_type { + auto delta = size() * b0_; + // Need to wrap in input space x, not output space y in case of a non-identify + // transform. + x -= std::floor((x - min_) / delta) * delta; + x += min_; + auto y = (this->forward(x) - min_) * b0_inv_; + return static_cast(y); + }; + + return detail::static_if>( + [&](auto x) { + if (std::isfinite(x)) + return lambda_circ(x); + else + return size(); + }, + lambda_circ, x); + } + index_type size_{0}; - internal_value_type min_{0}, delta_{1}; + internal_value_type min_{0}, b0_{0}, b0_inv_{0}; template friend class regular; From 4d9393e844f3258bc5e29900366d31c73e340b0f Mon Sep 17 00:00:00 2001 From: Ryan Date: Thu, 15 Jun 2023 19:00:20 -0400 Subject: [PATCH 2/5] added test for vakokako #336 example --- test/axis_regular_test.cpp | 38 ++++++++++++++++++++++++++++++++++++++ 1 file changed, 38 insertions(+) diff --git a/test/axis_regular_test.cpp b/test/axis_regular_test.cpp index 0e7cb0d6..67e682f6 100644 --- a/test/axis_regular_test.cpp +++ b/test/axis_regular_test.cpp @@ -301,5 +301,43 @@ int main() { BOOST_TEST_EQ(b.value(2), 5); } + // 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_regular = [](int N, double x0, double xN) { + return boost::histogram::axis::regular(N, x0, xN); + }; + + // The original example + fn_test_precision(27000, 0, 27000, fn_test_regular); + + // 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_regular); + } + } + } + } + return boost::report_errors(); } From e9b1a5ebeb1b36865fc8b8f52ab046d030d7fc07 Mon Sep 17 00:00:00 2001 From: Ryan Date: Thu, 15 Jun 2023 19:00:32 -0400 Subject: [PATCH 3/5] update serialization xml files --- test/axis_variant_serialization_test.xml | 5 +-- test/histogram_serialization_test_dynamic.xml | 31 +++++++++++-------- test/histogram_serialization_test_static.xml | 16 ++++++---- 3 files changed, 31 insertions(+), 21 deletions(-) diff --git a/test/axis_variant_serialization_test.xml b/test/axis_variant_serialization_test.xml index 9504c511..f17d4927 100644 --- a/test/axis_variant_serialization_test.xml +++ b/test/axis_variant_serialization_test.xml @@ -8,7 +8,7 @@ - + 1 @@ -17,7 +17,8 @@ 1 0.00000000000000000e+00 - 1.00000000000000000e+00 + 1.00000000000000000e+00 + 1.00000000000000000e+00 diff --git a/test/histogram_serialization_test_dynamic.xml b/test/histogram_serialization_test_dynamic.xml index 3411c758..26f5d04e 100644 --- a/test/histogram_serialization_test_dynamic.xml +++ b/test/histogram_serialization_test_dynamic.xml @@ -8,49 +8,52 @@ - + 7 0 - + 0 1 reg -1.00000000000000000e+00 - 2.00000000000000000e+00 + 2.00000000000000000e+00 + 5.00000000000000000e-01 - + 1 1 cir 0.000000000e+00 - 1.000000000e+00 + 1.000000000e+00 + 1.000000000e+00 - + 2 1 reg-log 0.00000000000000000e+00 - 2.00000000000000000e+00 + 2.88539008177792677e+00 + 3.46573590279972643e-01 - + 3 @@ -65,12 +68,13 @@ 3 1.00000000000000000e+00 - 9.00000000000000000e+00 + 9.00000000000000000e+00 + 1.11111111111111105e-01 - + 4 @@ -84,7 +88,7 @@ - + 5 @@ -98,7 +102,7 @@ - + 6 1 @@ -111,7 +115,7 @@ 0 4 - + 0 0 1 @@ -120,3 +124,4 @@ + diff --git a/test/histogram_serialization_test_static.xml b/test/histogram_serialization_test_static.xml index 95699780..6eb2216b 100644 --- a/test/histogram_serialization_test_static.xml +++ b/test/histogram_serialization_test_static.xml @@ -8,7 +8,7 @@ - + @@ -16,21 +16,24 @@ 1 reg -1.00000000000000000e+00 - 2.00000000000000000e+00 + 2.00000000000000000e+00 + 5.00000000000000000e-01 1 cir 0.000000000e+00 - 1.000000000e+00 + 1.000000000e+00 + 1.000000000e+00 1 reg-log 0.00000000000000000e+00 - 2.00000000000000000e+00 + 2.88539008177792677e+00 + 3.46573590279972643e-01 @@ -45,7 +48,8 @@ 3 1.00000000000000000e+00 - 9.00000000000000000e+00 + 9.00000000000000000e+00 + 1.11111111111111105e-01 @@ -74,7 +78,7 @@ 0 4 - + 0 0 1 From ab285420c509295c7c7254d26f8707183c4100d3 Mon Sep 17 00:00:00 2001 From: Ryan Date: Fri, 16 Jun 2023 10:19:00 -0400 Subject: [PATCH 4/5] fixed assert issue with -0.0 and +0.0 --- examples/getting_started_listing_01.cpp | 10 +++++++--- 1 file changed, 7 insertions(+), 3 deletions(-) diff --git a/examples/getting_started_listing_01.cpp b/examples/getting_started_listing_01.cpp index 7d74f380..66c027bb 100644 --- a/examples/getting_started_listing_01.cpp +++ b/examples/getting_started_listing_01.cpp @@ -72,18 +72,22 @@ int main() { Edges can be accessed with methods `lower()` and `upper()`. */ + // IEEE 754 has a positive and a negative zero, which are equal, but have different signs. + // This lambda function converts negative zero to positive zero for the purposes of string formatting. + auto fn_negative_zero_to_zero = [](double x) { return (x == 0) ? 0.0 : x; }; + std::ostringstream os; for (auto&& x : indexed(h, coverage::all)) { os << boost::format("bin %2i [%4.1f, %4.1f): %i\n") - % x.index() % x.bin().lower() % x.bin().upper() % *x; + % x.index() % fn_negative_zero_to_zero(x.bin().lower()) % fn_negative_zero_to_zero(x.bin().upper()) % *x; } std::cout << os.str() << std::flush; assert(os.str() == "bin -1 [-inf, -1.0): 1\n" "bin 0 [-1.0, -0.5): 1\n" - "bin 1 [-0.5, -0.0): 1\n" - "bin 2 [-0.0, 0.5): 2\n" + "bin 1 [-0.5, 0.0): 1\n" + "bin 2 [ 0.0, 0.5): 2\n" "bin 3 [ 0.5, 1.0): 0\n" "bin 4 [ 1.0, 1.5): 1\n" "bin 5 [ 1.5, 2.0): 1\n" From fb352576f20f0aa0b8d34d2a077966f3e2ca2654 Mon Sep 17 00:00:00 2001 From: Ryan Date: Fri, 16 Jun 2023 10:20:35 -0400 Subject: [PATCH 5/5] made lines at most 80 characters long --- examples/getting_started_listing_01.cpp | 8 +++++--- 1 file changed, 5 insertions(+), 3 deletions(-) diff --git a/examples/getting_started_listing_01.cpp b/examples/getting_started_listing_01.cpp index 66c027bb..146e442f 100644 --- a/examples/getting_started_listing_01.cpp +++ b/examples/getting_started_listing_01.cpp @@ -72,14 +72,16 @@ int main() { Edges can be accessed with methods `lower()` and `upper()`. */ - // IEEE 754 has a positive and a negative zero, which are equal, but have different signs. - // This lambda function converts negative zero to positive zero for the purposes of string formatting. + // IEEE 754 has a positive and a negative zero, which are equal, but have + // different signs. This lambda function converts negative zero to positive + // zero for the purposes of string formatting. auto fn_negative_zero_to_zero = [](double x) { return (x == 0) ? 0.0 : x; }; std::ostringstream os; for (auto&& x : indexed(h, coverage::all)) { os << boost::format("bin %2i [%4.1f, %4.1f): %i\n") - % x.index() % fn_negative_zero_to_zero(x.bin().lower()) % fn_negative_zero_to_zero(x.bin().upper()) % *x; + % x.index() % fn_negative_zero_to_zero(x.bin().lower()) + % fn_negative_zero_to_zero(x.bin().upper()) % *x; } std::cout << os.str() << std::flush;