Skip to content
New issue

Have a question about this project? Sign up for a free GitHub account to open an issue and contact its maintainers and the community.

Sign up for GitHub

By clicking “Sign up for GitHub”, you agree to our terms of service and privacy statement. We’ll occasionally send you account related emails.

Already on GitHub? Sign in to your account

Jump to bottom

Improve axis regular precision #387

Open
wants to merge 5 commits into
base: develop
Choose a base branch
from
Open

Improve axis regular precision #387

Show file tree
Hide file tree
Changes from all commits
Commits
Show all changes
5 commits
Select commit Hold shift + click to select a range
b2ac83c
refactored regular
ryanelandt Jun 15, 2023
4d9393e
added test for vakokako #336 example
ryanelandt Jun 15, 2023
e9b1a5e
update serialization xml files
ryanelandt Jun 15, 2023
ab28542
fixed assert issue with -0.0 and +0.0
ryanelandt Jun 16, 2023
fb35257
made lines at most 80 characters long
ryanelandt Jun 16, 2023
File filter

Filter by extension

Filter by extension
Conversations
Failed to load comments.
Loading
Jump to
Jump to file
Failed to load files.
Loading
Diff view
Diff view
12 changes: 9 additions & 3 deletions examples/getting_started_listing_01.cpp
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -72,18 +72,24 @@ 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"
Expand Down
120 changes: 71 additions & 49 deletions include/boost/histogram/axis/regular.hpp
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -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 <class T>
static T forward(T x) {
return std::log(x);
return std::log2(x);
}

/// Returns exp(x) for internal value x.
template <class T>
static T inverse(T x) {
return std::exp(x);
return std::pow(2, x);
}

template <class Archive>
Expand Down Expand Up @@ -221,8 +222,13 @@ class regular : public iterator_mixin<regular<Value, Transform, MetaData, Option
: transform_type(std::move(trans))
, metadata_base(std::move(meta))
, size_(static_cast<index_type>(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<internal_value_type>(delta / size());
b0_inv_ = static_cast<internal_value_type>(size() / delta);
}

// static_asserts were moved here from class scope to satisfy deduction in gcc>=11
static_assert(std::is_nothrow_move_constructible<transform_type>::value,
"transform must be no-throw move constructible");
Expand All @@ -233,11 +239,10 @@ class regular : public iterator_mixin<regular<Value, Transform, MetaData, Option
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_))
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).
Expand Down Expand Up @@ -306,52 +311,32 @@ class regular : public iterator_mixin<regular<Value, Transform, MetaData, Option

/// Return index for value argument.
index_type index(value_type x) const noexcept {
// Runs in hot loop, please measure impact of changes
auto z = (this->forward(x / unit_type{}) - min_) / delta_;
if (options_type::test(option::circular)) {
if (std::isfinite(z)) {
z -= std::floor(z);
return static_cast<index_type>(z * size());
}
} else {
if (z < 1) {
if (z >= 0)
return static_cast<index_type>(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<double>(x / unit_type{}));
}

/// Returns index and shift (if axis has grown) for the passed argument.
std::pair<index_type, index_type> 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<axis::index_type>(z * size());
auto y = (this->forward(x / unit_type{}) - min_) * b0_inv_;
if (y < size()) {
if (0 <= y) {
const auto i = static_cast<axis::index_type>(y);
return {i, 0};
}
if (z != -std::numeric_limits<internal_value_type>::infinity()) {
const auto stop = min_ + delta_;
const auto i = static_cast<axis::index_type>(std::floor(z * size()));
min_ += i * (delta_ / size());
delta_ = stop - min_;
if (y != -std::numeric_limits<internal_value_type>::infinity()) {
const auto i = static_cast<axis::index_type>(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<internal_value_type>::infinity()) {
const auto i = static_cast<axis::index_type>(z * size());
// y either beyond range, infinite, or NaN
if (y < std::numeric_limits<internal_value_type>::infinity()) {
const auto i = static_cast<axis::index_type>(y);
const auto n = i - size() + 1;
delta_ /= size();
delta_ *= size() + n;
size_ += n;
return {i, -n};
}
Expand All @@ -361,13 +346,13 @@ class regular : public iterator_mixin<regular<Value, Transform, MetaData, Option

/// Return value for fractional index argument.
value_type value(real_index_type i) const noexcept {
auto z = i / size();
if (!options_type::test(option::circular) && z < 0.0)
z = -std::numeric_limits<internal_value_type>::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<internal_value_type>::infinity() * b0_;
else if (options_type::test(option::circular) || i <= size())
z = min_ + i * b0_;
else {
z = std::numeric_limits<internal_value_type>::infinity() * delta_;
z = std::numeric_limits<internal_value_type>::infinity() * b0_;
}
return static_cast<value_type>(this->inverse(z) * unit_type());
}
Expand All @@ -386,7 +371,7 @@ class regular : public iterator_mixin<regular<Value, Transform, MetaData, Option
template <class V, class T, class M, class O>
bool operator==(const regular<V, T, M, O>& 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 <class V, class T, class M, class O>
Expand All @@ -400,12 +385,49 @@ class regular : public iterator_mixin<regular<Value, Transform, MetaData, Option
ar& make_nvp("size", size_);
ar& make_nvp("meta", this->metadata());
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<index_type>(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<index_type>(y);
};

return detail::static_if<std::is_floating_point<value_type>>(
[&](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 <class V, class T, class M, class O>
friend class regular;
Expand Down
38 changes: 38 additions & 0 deletions test/axis_regular_test.cpp
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -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<double>(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<double> v_spacing = {0.125, 0.375, 0.5, 0.75, 1, 1.625, 3, 7.25};
const std::vector<double> 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();
}
5 changes: 3 additions & 2 deletions test/axis_variant_serialization_test.xml
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -8,7 +8,7 @@

<?xml version="1.0" encoding="UTF-8" standalone="yes" ?>
<!DOCTYPE boost_serialization>
<boost_serialization signature="serialization::archive" version="17">
<boost_serialization signature="serialization::archive" version="19">
<item class_id="0" tracking_level="0" version="0">
<variant class_id="1" tracking_level="0" version="0">
<which>1</which>
Expand All @@ -17,7 +17,8 @@
<size>1</size>
<meta></meta>
<min>0.00000000000000000e+00</min>
<delta>1.00000000000000000e+00</delta>
<b0>1.00000000000000000e+00</b0>
<b0_inv>1.00000000000000000e+00</b0_inv>
</value>
</variant>
</item>
Expand Down
31 changes: 18 additions & 13 deletions test/histogram_serialization_test_dynamic.xml
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -8,49 +8,52 @@

<?xml version="1.0" encoding="UTF-8" standalone="yes" ?>
<!DOCTYPE boost_serialization>
<boost_serialization signature="serialization::archive" version="17">
<boost_serialization signature="serialization::archive" version="19">
<item class_id="0" tracking_level="0" version="0">
<axes class_id="1" tracking_level="0" version="0">
<count>7</count>
<item_version>0</item_version>
<item class_id="2" tracking_level="0" version="0">
<variant class_id="3" tracking_level="1" version="0" object_id="_0">
<variant class_id="3" tracking_level="0" version="0">
<which>0</which>
<value class_id="4" tracking_level="0" version="0">
<transform class_id="5" tracking_level="0" version="0"></transform>
<size>1</size>
<meta>reg</meta>
<min>-1.00000000000000000e+00</min>
<delta>2.00000000000000000e+00</delta>
<b0>2.00000000000000000e+00</b0>
<b0_inv>5.00000000000000000e-01</b0_inv>
</value>
</variant>
</item>
<item>
<variant object_id="_1">
<variant>
<which>1</which>
<value class_id="6" tracking_level="0" version="0">
<transform></transform>
<size>1</size>
<meta>cir</meta>
<min>0.000000000e+00</min>
<delta>1.000000000e+00</delta>
<b0>1.000000000e+00</b0>
<b0_inv>1.000000000e+00</b0_inv>
</value>
</variant>
</item>
<item>
<variant object_id="_2">
<variant>
<which>2</which>
<value class_id="7" tracking_level="0" version="0">
<transform class_id="8" tracking_level="0" version="0"></transform>
<size>1</size>
<meta>reg-log</meta>
<min>0.00000000000000000e+00</min>
<delta>2.00000000000000000e+00</delta>
<b0>2.88539008177792677e+00</b0>
<b0_inv>3.46573590279972643e-01</b0_inv>
</value>
</variant>
</item>
<item>
<variant object_id="_3">
<variant>
<which>3</which>
<value class_id="9" tracking_level="0" version="0">
<transform class_id="10" tracking_level="0" version="0">
Expand All @@ -65,12 +68,13 @@
<item>3</item>
</meta>
<min>1.00000000000000000e+00</min>
<delta>9.00000000000000000e+00</delta>
<b0>9.00000000000000000e+00</b0>
<b0_inv>1.11111111111111105e-01</b0_inv>
</value>
</variant>
</item>
<item>
<variant object_id="_4">
<variant>
<which>4</which>
<value class_id="12" tracking_level="0" version="0">
<seq>
Expand All @@ -84,7 +88,7 @@
</variant>
</item>
<item>
<variant object_id="_5">
<variant>
<which>5</which>
<value class_id="14" tracking_level="0" version="0">
<seq>
Expand All @@ -98,7 +102,7 @@
</variant>
</item>
<item>
<variant object_id="_6">
<variant>
<which>6</which>
<value class_id="15" tracking_level="0" version="0">
<size>1</size>
Expand All @@ -111,7 +115,7 @@
<storage class_id="17" tracking_level="0" version="0">
<type>0</type>
<size>4</size>
<buffer>
<buffer class_id="18" tracking_level="0" version="0">
<item>0</item>
<item>0</item>
<item>1</item>
Expand All @@ -120,3 +124,4 @@
</storage>
</item>
</boost_serialization>

Loading