Fix scalar-broadcast fill after a non-inclusive axis (scikit-hep/boost-histogram#960)

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 <class T>
-  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<T, value_type>::value ||
+                          !is_iterable<T>::value)>{},
+           iterable_or_value);
+  }
+
+  template <class T>
+  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 <class T>
+  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,43 +101,60 @@ struct index_visitor {
   }
 
   template <class T>
-  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<value_type, std::invalid_argument>(x));
   }
 
   template <class T>
-  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 <class T>
-  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
 
-    // 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<std::intptr_t>(*begin_) - static_cast<std::intptr_t>(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_};
+    linearize_scalar(IsGrowing{}, idx, value);
+    if (is_valid(idx)) {
+      // 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<std::intptr_t>(static_cast<std::size_t>(idx)) -
+                         static_cast<std::intptr_t>(stride_);
+      for (auto it = begin_; it != begin_ + size_; ++it) *it += delta;
+    } else {
       std::fill(begin_, begin_ + size_, invalid_index);
+    }
   }
 
   template <class T>
-  void operator()(const T& iterable_or_value) const {
-    call_1(mp11::mp_bool<(std::is_convertible<T, value_type>::value ||
-                          !is_iterable<T>::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_type, std::invalid_argument>(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 <class T>
+  void linearize_scalar(std::false_type, index_type& idx, const T& value) const {
+    // no axis is growing
+    linearize(idx, stride_, axis_, try_cast<value_type, std::invalid_argument>(value));
   }
+
 };
 
 template <class Index, class S, class Axes, class T>

test/histogram_fill_test.cpp

@@ -44,7 +44,10 @@ using ing = axis::integer<double, axis::null_type,
                           decltype(axis::option::growth | axis::option::underflow |
                                    axis::option::overflow)>;
 using cs = axis::category<std::string, axis::null_type>;
+using cs0 = axis::category<std::string, axis::null_type, axis::option::none_t>;
 using csg = axis::category<std::string, axis::null_type, axis::option::growth_t>;
+using ci = axis::category<int, axis::null_type>;
+using ci0 = axis::category<int, axis::null_type, axis::option::none_t>;
 
 struct axis2d {
   auto size() const { return axis::index_type{2}; }
@@ -345,6 +348,191 @@ void run_tests(const std::vector<int>& x, const std::vector<int>& 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 <class Tag>
+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<int> 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<int, std::vector<int>, std::string, std::vector<std::string>>;
+    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<int> data = {4, 1, 2, 3};
+    std::vector<double> 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<int, std::vector<int>, std::string, std::vector<std::string>>;
+    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<int> data = {4, 1, 2, 3};
+
+    for (auto&& xi : data) h(xi, "A");
+
+    using V = variant<int, std::vector<int>, std::string, std::vector<std::string>>;
+    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<int> data1 = {99, 1, 2};
+    for (auto&& xi : data1) h(xi, "A");
+    using V = variant<int, std::vector<int>, 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<int> 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<int> 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<std::string> data = {"C", "A", "B"}; // "C" is out-of-bounds
+
+    for (auto&& si : data) h(1, si);
+
+    using V = variant<int, std::vector<int>, std::string, std::vector<std::string>>;
+    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<int> data = {0, 1, 2}; // 0 is out-of-bounds for in0{1, 3}
+
+    for (auto&& xi : data) h(xi, "A");
+
+    using V = variant<int, std::vector<int>, std::string, std::vector<std::string>>;
+    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<int, axis::null_type, axis::option::growth_t>;
+    auto h = make(Tag(), in0{1, 3}, ig{0, 2});
+    auto h2 = h;
+
+    std::vector<int> 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<int, std::vector<int>>;
+    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<int> 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<int, std::vector<int>>;
+    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<static_tag>(x, y, w);
   run_tests<dynamic_tag>(x, y, w);
 
+  run_tests_issue_960<static_tag>();
+  run_tests_issue_960<dynamic_tag>();
+
   return boost::report_errors();
 }