From 8232088c159d8c38c1218a9000f160327cf866fc Mon Sep 17 00:00:00 2001 From: Nick Manganelli Date: Tue, 14 Jul 2026 07:43:25 +0200 Subject: [PATCH 1/6] Fix scalar-broadcast fill after a non-inclusive axis (#960) In the vectorized fill (fill_n), a scalar argument is broadcast across an axis by computing that axis' contribution to the linear index once and adding it to every entry. The contribution was derived from the change of the first index (*begin_). When a previous axis had already invalidated some entries, including the first one, the broadcast axis then invalidated *all* entries, silently discarding valid data. This is scikit-hep/boost-histogram#960: a non-growing axis combined with out-of-bounds entries and a broadcast scalar fill produced a sum of zero. Compute the axis contribution on a fresh index instead and add it to all entries, leaving already-invalid entries untouched. The growing and non-growing scalar paths share the same code via linearize/linearize_growth. Co-Authored-By: Claude Opus 4.8 --- include/boost/histogram/detail/fill_n.hpp | 42 ++++++++++++++++------- 1 file changed, 30 insertions(+), 12 deletions(-) diff --git a/include/boost/histogram/detail/fill_n.hpp b/include/boost/histogram/detail/fill_n.hpp index e5a540bc..ec73ff59 100644 --- a/include/boost/histogram/detail/fill_n.hpp +++ b/include/boost/histogram/detail/fill_n.hpp @@ -103,19 +103,37 @@ struct index_visitor { void call_1(std::true_type, const T& value) const { // T is compatible value; fill single value N times - // Optimization: We call call_2 only once and then add the index shift onto the - // whole array of indices, because it is always the same. This also works if the - // axis grows during this operation. There are no shifts to apply if the zero-point - // changes. - const auto before = *begin_; - call_2(IsGrowing{}, begin_, value); - if (is_valid(*begin_)) { - // since index can be std::size_t or optional_index, must do conversion here - const auto delta = - static_cast(*begin_) - static_cast(before); - for (auto it = begin_ + 1; it != begin_ + size_; ++it) *it += delta; - } else + // Compute the contribution of this axis to the linear index once on a + // fresh index and add it to the whole array of indices. The contribution + // must not be derived from the change of *begin_, because a previous axis + // may have already invalidated some entries, including the first one + // (scikit-hep/boost-histogram#960). Entries that are already invalid stay + // invalid, all other entries receive the same contribution. + // + // The fresh index is biased by stride_, because the contribution of an + // underflow bin is -stride_ and the index type is unsigned. The bias is + // subtracted again when the delta is computed. + index_type idx{stride_}; + if (IsGrowing::value) { + // IsGrowing::value is a compile-time constant, the dead branch is + // eliminated; linearize_growth also handles non-growing axes. + axis::index_type shift; + linearize_growth(idx, shift, stride_, axis_, + try_cast(value)); + // No index shifts to apply if the zero-point changes, since all entries + // receive the same contribution which is computed after the growth, but + // the shift must be recorded so that the storage is resized correctly. + if (shift > 0) *shift_ += shift; + } else { + linearize(idx, stride_, axis_, try_cast(value)); + } + if (is_valid(idx)) { + const auto delta = static_cast(static_cast(idx)) - + static_cast(stride_); + for (auto it = begin_; it != begin_ + size_; ++it) *it += delta; + } else { std::fill(begin_, begin_ + size_, invalid_index); + } } template From 9695104442e4f1c9d883b1a5684594a8833ddc42 Mon Sep 17 00:00:00 2001 From: Nick Manganelli Date: Tue, 14 Jul 2026 07:43:25 +0200 Subject: [PATCH 2/6] Add regression tests for scikit-hep/boost-histogram#960 Compare the vectorized fill against the one-at-a-time reference for scalar broadcast after a non-inclusive axis, covering the growing-axis, zero-point-shift, and underflow variants. Co-Authored-By: Claude Opus 4.8 --- test/histogram_fill_test.cpp | 191 +++++++++++++++++++++++++++++++++++ 1 file changed, 191 insertions(+) diff --git a/test/histogram_fill_test.cpp b/test/histogram_fill_test.cpp index d647412f..9d97f7a0 100644 --- a/test/histogram_fill_test.cpp +++ b/test/histogram_fill_test.cpp @@ -44,7 +44,10 @@ using ing = axis::integer; using cs = axis::category; +using cs0 = axis::category; using csg = axis::category; +using ci = axis::category; +using ci0 = axis::category; struct axis2d { auto size() const { return axis::index_type{2}; } @@ -345,6 +348,191 @@ void run_tests(const std::vector& x, const std::vector& y, } } +// Test for GitHub issue scikit-hep/boost-histogram#960: +// Filling with scalar value for one axis and array for another, where array contains +// out-of-bounds values and axis has no overflow, should correctly accumulate valid +// entries instead of producing sum of zero. +template +void run_tests_issue_960() { + // 2D histogram with IntCategory (no overflow) and StrCategory + // Bug: when filling with array for first axis (containing out-of-bounds values) + // and scalar for second axis, valid entries were incorrectly discarded. + { + auto h = make(Tag(), ci0{1, 2, 3}, cs{"A", "B"}); + auto h2 = h; + + // Data with out-of-bounds value (4) and valid values (1, 2, 3) + std::vector data = {4, 1, 2, 3}; + + // Fill one-at-a-time (reference) + for (auto&& xi : data) h(xi, "A"); + + // Fill using array for first axis, scalar for second axis + using V = variant, std::string, std::vector>; + V xy[2]; + xy[0] = data; + xy[1] = "A"; + h2.fill(xy); + + BOOST_TEST_EQ(h, h2); + BOOST_TEST_EQ(sum(h), 3); // only 1, 2, 3 should be counted + } + + // Same test with weight storage + { + auto h = make_s(Tag(), weight_storage(), ci0{1, 2, 3}, cs{"A", "B"}); + auto h2 = h; + + std::vector data = {4, 1, 2, 3}; + std::vector w = {1.0, 1.0, 1.0, 1.0}; + + for (unsigned i = 0; i < data.size(); ++i) h(data[i], "A", weight(w[i])); + + using V = variant, std::string, std::vector>; + V xy[2]; + xy[0] = data; + xy[1] = "A"; + h2.fill(xy, weight(w)); + + BOOST_TEST_EQ(h, h2); + BOOST_TEST_EQ(sum(h).value(), 3); // only 1, 2, 3 should be counted + } + + // Test with both axes having no overflow + { + auto h = make(Tag(), ci0{1, 2, 3}, cs0{"A", "B"}); + auto h2 = h; + + std::vector data = {4, 1, 2, 3}; + + for (auto&& xi : data) h(xi, "A"); + + using V = variant, std::string, std::vector>; + V xy[2]; + xy[0] = data; + xy[1] = "A"; + h2.fill(xy); + + BOOST_TEST_EQ(h, h2); + BOOST_TEST_EQ(sum(h), 3); + } + + // Test with out-of-bounds value at different positions in array + { + auto h = make(Tag(), ci0{1, 2, 3}, cs{"A"}); + auto h2 = h; + + // Out-of-bounds at start + std::vector data1 = {99, 1, 2}; + for (auto&& xi : data1) h(xi, "A"); + using V = variant, std::string>; + V xy[2]; + xy[0] = data1; + xy[1] = "A"; + h2.fill(xy); + BOOST_TEST_EQ(h, h2); + + // Out-of-bounds in middle + std::vector data2 = {1, 99, 2}; + for (auto&& xi : data2) h(xi, "A"); + xy[0] = data2; + h2.fill(xy); + BOOST_TEST_EQ(h, h2); + + // Out-of-bounds at end + std::vector data3 = {1, 2, 99}; + for (auto&& xi : data3) h(xi, "A"); + xy[0] = data3; + h2.fill(xy); + BOOST_TEST_EQ(h, h2); + + BOOST_TEST_EQ(sum(h), 6); // 2 valid values * 3 fills + } + + // Test scalar for first axis, array for second axis + { + auto h = make(Tag(), ci0{1, 2, 3}, cs0{"A", "B"}); + auto h2 = h; + + std::vector data = {"C", "A", "B"}; // "C" is out-of-bounds + + for (auto&& si : data) h(1, si); + + using V = variant, std::string, std::vector>; + V xy[2]; + xy[0] = 1; + xy[1] = data; + h2.fill(xy); + + BOOST_TEST_EQ(h, h2); + BOOST_TEST_EQ(sum(h), 2); // only "A" and "B" should be counted + } + + // Scalar broadcast on a growing axis, while the previous non-inclusive axis + // has already invalidated the first entry (growing variant of the bug) + { + auto h = make(Tag(), in0{1, 3}, csg{"X", "A"}); + auto h2 = h; + + std::vector data = {0, 1, 2}; // 0 is out-of-bounds for in0{1, 3} + + for (auto&& xi : data) h(xi, "A"); + + using V = variant, std::string, std::vector>; + V xy[2]; + xy[0] = data; + xy[1] = "A"; + h2.fill(xy); + + BOOST_TEST_EQ(h, h2); + BOOST_TEST_EQ(sum(h), 2); + } + + // Scalar broadcast that triggers growth at the lower end (zero-point shift), + // while the previous non-inclusive axis has already invalidated the first entry + { + using ig = axis::integer; + auto h = make(Tag(), in0{1, 3}, ig{0, 2}); + auto h2 = h; + + std::vector data = {0, 1, 2}; // 0 is out-of-bounds for in0{1, 3} + + for (auto&& xi : data) h(xi, -1); // -1 grows ig downward + for (auto&& xi : data) h(xi, 5); // 5 grows ig upward + + using V = variant>; + V xy[2]; + xy[0] = data; + xy[1] = -1; + h2.fill(xy); + xy[1] = 5; + h2.fill(xy); + + BOOST_TEST_EQ(h, h2); + BOOST_TEST_EQ(sum(h), 4); + } + + // Scalar broadcast that lands in the underflow bin (negative index + // contribution), while the first entry is already invalid + { + auto h = make(Tag(), in0{1, 3}, in{1, 3}); + auto h2 = h; + + std::vector data = {0, 1, 2}; // 0 is out-of-bounds for in0{1, 3} + + for (auto&& xi : data) h(xi, 0); // 0 lands in underflow of in{1, 3} + + using V = variant>; + V xy[2]; + xy[0] = data; + xy[1] = 0; + h2.fill(xy); + + BOOST_TEST_EQ(h, h2); + BOOST_TEST_EQ(sum(h), 2); + } +} + int main() { std::mt19937 gen(1); std::normal_distribution<> id(0, 2); @@ -359,5 +547,8 @@ int main() { run_tests(x, y, w); run_tests(x, y, w); + run_tests_issue_960(); + run_tests_issue_960(); + return boost::report_errors(); } From 98d428d83ad66e85435944b9a7c0ea32f0228a54 Mon Sep 17 00:00:00 2001 From: Hans Dembinski Date: Sun, 11 Oct 2026 11:23:03 +0000 Subject: [PATCH 3/6] Address review comments on #431 - Remove call_2, inlining the per-element linearize calls into the iterable branch of the visitor - Rename call_1 to call_impl - Explain why the double cast in the scalar-broadcast delta is intentional, not superfluous --- include/boost/histogram/detail/fill_n.hpp | 49 +++++++++++++---------- 1 file changed, 27 insertions(+), 22 deletions(-) diff --git a/include/boost/histogram/detail/fill_n.hpp b/include/boost/histogram/detail/fill_n.hpp index ec73ff59..ece5dc76 100644 --- a/include/boost/histogram/detail/fill_n.hpp +++ b/include/boost/histogram/detail/fill_n.hpp @@ -75,32 +75,30 @@ struct index_visitor { : axis_(a), stride_(str), start_(sta), size_(si), begin_(it), shift_(shift) {} template - void call_2(std::true_type, pointer it, const T& x) const { - // must use this code for all axes if one of them is growing - axis::index_type shift; - linearize_growth(*it, shift, stride_, axis_, - try_cast(x)); - if (shift > 0) { // shift previous indices, because axis zero-point has changed - while (it != begin_) *--it += static_cast(shift) * stride_; - *shift_ += shift; - } - } - - template - void call_2(std::false_type, pointer it, const T& x) const { - // no axis is growing - linearize(*it, stride_, axis_, try_cast(x)); - } - - template - void call_1(std::false_type, const T& iterable) const { + void call_impl(std::false_type, const T& iterable) const { // T is iterable; fill N values const auto* tp = dtl::data(iterable) + start_; - for (auto it = begin_; it != begin_ + size_; ++it) call_2(IsGrowing{}, it, *tp++); + if (IsGrowing::value) { + // must use this code for all axes if one of them is growing + for (auto it = begin_; it != begin_ + size_; ++it) { + axis::index_type shift; + linearize_growth(*it, shift, stride_, axis_, + try_cast(*tp++)); + if (shift > 0) { // shift previous indices, because axis zero-point has changed + while (it != begin_) *--it += static_cast(shift) * stride_; + *shift_ += shift; + } + } + } else { + // no axis is growing + for (auto it = begin_; it != begin_ + size_; ++it) + linearize(*it, stride_, axis_, + try_cast(*tp++)); + } } template - void call_1(std::true_type, const T& value) const { + void call_impl(std::true_type, const T& value) const { // T is compatible value; fill single value N times // Compute the contribution of this axis to the linear index once on a @@ -128,6 +126,13 @@ struct index_visitor { linearize(idx, stride_, axis_, try_cast(value)); } if (is_valid(idx)) { + // The double cast is intentional, not superfluous. The inner std::size_t + // cast extracts the raw value, because Index is optional_index for + // growing axes and has no direct cast to std::intptr_t. The outer cast + // makes the subtraction signed, so that the underflow contribution + // (-stride_) yields a negative delta instead of wrapping to a huge + // positive value. Both values are bounded by the storage size, so the + // conversion is value-preserving, not a two's complement trick. const auto delta = static_cast(static_cast(idx)) - static_cast(stride_); for (auto it = begin_; it != begin_ + size_; ++it) *it += delta; @@ -138,7 +143,7 @@ struct index_visitor { template void operator()(const T& iterable_or_value) const { - call_1(mp11::mp_bool<(std::is_convertible::value || + call_impl(mp11::mp_bool<(std::is_convertible::value || !is_iterable::value)>{}, iterable_or_value); } From cb26ca3e6e0e7ab68879d243cce7aad48bce891e Mon Sep 17 00:00:00 2001 From: Hans Dembinski Date: Sun, 11 Oct 2026 14:13:09 +0200 Subject: [PATCH 4/6] Revert rename of call_1 and inlining of call_2, keep double-cast comment --- include/boost/histogram/detail/fill_n.hpp | 42 ++++++++++++----------- 1 file changed, 22 insertions(+), 20 deletions(-) diff --git a/include/boost/histogram/detail/fill_n.hpp b/include/boost/histogram/detail/fill_n.hpp index ece5dc76..6666e3f9 100644 --- a/include/boost/histogram/detail/fill_n.hpp +++ b/include/boost/histogram/detail/fill_n.hpp @@ -75,30 +75,32 @@ struct index_visitor { : axis_(a), stride_(str), start_(sta), size_(si), begin_(it), shift_(shift) {} template - void call_impl(std::false_type, const T& iterable) const { + void call_2(std::true_type, pointer it, const T& x) const { + // must use this code for all axes if one of them is growing + axis::index_type shift; + linearize_growth(*it, shift, stride_, axis_, + try_cast(x)); + if (shift > 0) { // shift previous indices, because axis zero-point has changed + while (it != begin_) *--it += static_cast(shift) * stride_; + *shift_ += shift; + } + } + + template + void call_2(std::false_type, pointer it, const T& x) const { + // no axis is growing + linearize(*it, stride_, axis_, try_cast(x)); + } + + template + void call_1(std::false_type, const T& iterable) const { // T is iterable; fill N values const auto* tp = dtl::data(iterable) + start_; - if (IsGrowing::value) { - // must use this code for all axes if one of them is growing - for (auto it = begin_; it != begin_ + size_; ++it) { - axis::index_type shift; - linearize_growth(*it, shift, stride_, axis_, - try_cast(*tp++)); - if (shift > 0) { // shift previous indices, because axis zero-point has changed - while (it != begin_) *--it += static_cast(shift) * stride_; - *shift_ += shift; - } - } - } else { - // no axis is growing - for (auto it = begin_; it != begin_ + size_; ++it) - linearize(*it, stride_, axis_, - try_cast(*tp++)); - } + for (auto it = begin_; it != begin_ + size_; ++it) call_2(IsGrowing{}, it, *tp++); } template - void call_impl(std::true_type, const T& value) const { + void call_1(std::true_type, const T& value) const { // T is compatible value; fill single value N times // Compute the contribution of this axis to the linear index once on a @@ -143,7 +145,7 @@ struct index_visitor { template void operator()(const T& iterable_or_value) const { - call_impl(mp11::mp_bool<(std::is_convertible::value || + call_1(mp11::mp_bool<(std::is_convertible::value || !is_iterable::value)>{}, iterable_or_value); } From 575d396a07da0a2ec206d189798505e19868b13b Mon Sep 17 00:00:00 2001 From: Hans Dembinski Date: Sun, 11 Oct 2026 14:17:54 +0200 Subject: [PATCH 5/6] Restore static dispatch for scalar linearize in fill_n --- include/boost/histogram/detail/fill_n.hpp | 32 ++++++++++++++--------- 1 file changed, 19 insertions(+), 13 deletions(-) diff --git a/include/boost/histogram/detail/fill_n.hpp b/include/boost/histogram/detail/fill_n.hpp index 6666e3f9..4b4fc1cd 100644 --- a/include/boost/histogram/detail/fill_n.hpp +++ b/include/boost/histogram/detail/fill_n.hpp @@ -99,6 +99,24 @@ struct index_visitor { for (auto it = begin_; it != begin_ + size_; ++it) call_2(IsGrowing{}, it, *tp++); } + template + void scalar_linearize(std::true_type, index_type& idx, const T& value) const { + // must use this code for all axes if one of them is growing + axis::index_type shift; + linearize_growth(idx, shift, stride_, axis_, + try_cast(value)); + // No index shifts to apply if the zero-point changes, since all entries + // receive the same contribution which is computed after the growth, but + // the shift must be recorded so that the storage is resized correctly. + if (shift > 0) *shift_ += shift; + } + + template + void scalar_linearize(std::false_type, index_type& idx, const T& value) const { + // no axis is growing + linearize(idx, stride_, axis_, try_cast(value)); + } + template void call_1(std::true_type, const T& value) const { // T is compatible value; fill single value N times @@ -114,19 +132,7 @@ struct index_visitor { // underflow bin is -stride_ and the index type is unsigned. The bias is // subtracted again when the delta is computed. index_type idx{stride_}; - if (IsGrowing::value) { - // IsGrowing::value is a compile-time constant, the dead branch is - // eliminated; linearize_growth also handles non-growing axes. - axis::index_type shift; - linearize_growth(idx, shift, stride_, axis_, - try_cast(value)); - // No index shifts to apply if the zero-point changes, since all entries - // receive the same contribution which is computed after the growth, but - // the shift must be recorded so that the storage is resized correctly. - if (shift > 0) *shift_ += shift; - } else { - linearize(idx, stride_, axis_, try_cast(value)); - } + scalar_linearize(IsGrowing{}, idx, value); if (is_valid(idx)) { // The double cast is intentional, not superfluous. The inner std::size_t // cast extracts the raw value, because Index is optional_index for From b0b0cd36c62cc92815f42d4404ed32aa81e652df Mon Sep 17 00:00:00 2001 From: Hans Dembinski Date: Sun, 11 Oct 2026 14:34:25 +0200 Subject: [PATCH 6/6] Rename fill_n index_visitor helpers for clarity - call_1 -> dispatch_on_value_or_iterable - call_2 -> linearize_element (per-element linearize, pointer based) - scalar_linearize -> linearize_scalar (scalar broadcast, index by reference) - move operator() to the top and reorder the helpers --- include/boost/histogram/detail/fill_n.hpp | 80 +++++++++++------------ 1 file changed, 40 insertions(+), 40 deletions(-) diff --git a/include/boost/histogram/detail/fill_n.hpp b/include/boost/histogram/detail/fill_n.hpp index 4b4fc1cd..031b3d46 100644 --- a/include/boost/histogram/detail/fill_n.hpp +++ b/include/boost/histogram/detail/fill_n.hpp @@ -75,7 +75,21 @@ struct index_visitor { : axis_(a), stride_(str), start_(sta), size_(si), begin_(it), shift_(shift) {} template - void call_2(std::true_type, pointer it, const T& x) const { + void operator()(const T& iterable_or_value) const { + dispatch_on_value_or_iterable(mp11::mp_bool<(std::is_convertible::value || + !is_iterable::value)>{}, + iterable_or_value); + } + + template + void dispatch_on_value_or_iterable(std::false_type, const T& iterable) const { + // T is iterable; fill N values + const auto* tp = dtl::data(iterable) + start_; + for (auto it = begin_; it != begin_ + size_; ++it) linearize_element(IsGrowing{}, it, *tp++); + } + + template + void linearize_element(std::true_type, pointer it, const T& x) const { // must use this code for all axes if one of them is growing axis::index_type shift; linearize_growth(*it, shift, stride_, axis_, @@ -87,38 +101,13 @@ struct index_visitor { } template - void call_2(std::false_type, pointer it, const T& x) const { + void linearize_element(std::false_type, pointer it, const T& x) const { // no axis is growing linearize(*it, stride_, axis_, try_cast(x)); } template - void call_1(std::false_type, const T& iterable) const { - // T is iterable; fill N values - const auto* tp = dtl::data(iterable) + start_; - for (auto it = begin_; it != begin_ + size_; ++it) call_2(IsGrowing{}, it, *tp++); - } - - template - void scalar_linearize(std::true_type, index_type& idx, const T& value) const { - // must use this code for all axes if one of them is growing - axis::index_type shift; - linearize_growth(idx, shift, stride_, axis_, - try_cast(value)); - // No index shifts to apply if the zero-point changes, since all entries - // receive the same contribution which is computed after the growth, but - // the shift must be recorded so that the storage is resized correctly. - if (shift > 0) *shift_ += shift; - } - - template - void scalar_linearize(std::false_type, index_type& idx, const T& value) const { - // no axis is growing - linearize(idx, stride_, axis_, try_cast(value)); - } - - template - void call_1(std::true_type, const T& value) const { + void dispatch_on_value_or_iterable(std::true_type, const T& value) const { // T is compatible value; fill single value N times // Compute the contribution of this axis to the linear index once on a @@ -132,15 +121,14 @@ struct index_visitor { // underflow bin is -stride_ and the index type is unsigned. The bias is // subtracted again when the delta is computed. index_type idx{stride_}; - scalar_linearize(IsGrowing{}, idx, value); + linearize_scalar(IsGrowing{}, idx, value); if (is_valid(idx)) { - // The double cast is intentional, not superfluous. The inner std::size_t - // cast extracts the raw value, because Index is optional_index for - // growing axes and has no direct cast to std::intptr_t. The outer cast - // makes the subtraction signed, so that the underflow contribution - // (-stride_) yields a negative delta instead of wrapping to a huge - // positive value. Both values are bounded by the storage size, so the - // conversion is value-preserving, not a two's complement trick. + // The double cast is intentional. The inner std::size_t cast extracts the raw + // value, because Index is optional_index for growing axes and has no direct cast to + // std::intptr_t. The outer cast makes the subtraction signed, so that the underflow + // contribution (-stride_) yields a negative delta instead of wrapping to a huge + // positive value. Both values are bounded by the storage size, so the conversion is + // value-preserving, not a two's complement trick. const auto delta = static_cast(static_cast(idx)) - static_cast(stride_); for (auto it = begin_; it != begin_ + size_; ++it) *it += delta; @@ -150,11 +138,23 @@ struct index_visitor { } template - void operator()(const T& iterable_or_value) const { - call_1(mp11::mp_bool<(std::is_convertible::value || - !is_iterable::value)>{}, - iterable_or_value); + void linearize_scalar(std::true_type, index_type& idx, const T& value) const { + // must use this code for all axes if one of them is growing + axis::index_type shift; + linearize_growth(idx, shift, stride_, axis_, + try_cast(value)); + // No index shifts to apply if the zero-point changes, since all entries + // receive the same contribution which is computed after the growth, but + // the shift must be recorded so that the storage is resized correctly. + if (shift > 0) *shift_ += shift; } + + template + void linearize_scalar(std::false_type, index_type& idx, const T& value) const { + // no axis is growing + linearize(idx, stride_, axis_, try_cast(value)); + } + }; template