Fix bug in 1D histogram::fill when axis is growing and argument is single value (#334)

.github/workflows/slow.yml

@@ -64,9 +64,12 @@ jobs:
         ../../b2 $B2_OPTS cxxstd=17 test//all
 
   gcc5:
-    runs-on: ubuntu-16.04
+    runs-on: ubuntu-18.04
     steps:
     - uses: actions/checkout@v2
+    - uses: egor-tensin/setup-gcc@v1
+      with:
+        version: 5
     - name: Fetch Boost superproject
       run: |
         cd ..

include/boost/histogram/detail/axes.hpp

@@ -349,7 +349,21 @@ std::size_t bincount(const T& axes) {
   return n;
 }
 
-// initial offset for the linear index
+/** Initial offset for the linear index.
+
+  This precomputes an offset for the global index so that axis index = -1 addresses the
+  first entry in the storage. Example: one-dim. histogram. The offset is 1, stride is 1,
+  and global_index = offset + axis_index * stride == 0 addresses the first element of
+  the storage.
+
+  Using the offset makes the hot inner loop that computes the global_index simpler and
+  thus faster, because we do not have to branch for each axis to check whether it has
+  an underflow bin.
+
+  The offset is set to an invalid value when the histogram contains at least one growing
+  axis, because this optimization then cannot be used. See detail/linearize.hpp, in this
+  case linearize_growth is called.
+*/
 template <class T>
 std::size_t offset(const T& axes) {
   std::size_t n = 0;

include/boost/histogram/detail/fill.hpp

@@ -250,7 +250,7 @@ template <class ArgTraits, class Storage, class Axes, class Args>
 auto fill_2(ArgTraits, mp11::mp_true, const std::size_t, Storage& st, Axes& axes,
             const Args& args) {
   std::array<axis::index_type, ArgTraits::nargs::value> shifts;
-  // offset must be zero for linearize_growth
+  // offset must be zero for linearize_growth (value of offset argument is ignored)
   mp11::mp_if<has_non_inclusive_axis<Axes>, optional_index, std::size_t> idx{0};
   std::size_t stride = 1;
   bool update_needed = false;

include/boost/histogram/detail/fill_n.hpp

@@ -101,12 +101,18 @@ struct index_visitor {
   template <class T>
   void call_1(std::true_type, const T& value) const {
     // T is compatible value; fill single value N times
-    index_type idx{*begin_};
+
-    call_2(IsGrowing{}, &idx, value);
+    // Optimization: We call call_2 only once and then add the index shift onto the
-    if (is_valid(idx)) {
+    // 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>(idx) - static_cast<std::intptr_t>(*begin_);
+          static_cast<std::intptr_t>(*begin_) - static_cast<std::intptr_t>(before);
-      for (auto&& i : make_span(begin_, size_)) i += delta;
+      for (auto it = begin_ + 1; it != begin_ + size_; ++it) *it += delta;
     } else
       std::fill(begin_, begin_ + size_, invalid_index);
   }
@@ -129,7 +135,9 @@ void fill_n_indices(Index* indices, const std::size_t start, const std::size_t s
     *eit++ = axis::traits::extent(a);
   }); // LCOV_EXCL_LINE: gcc-8 is missing this line for no reason
 
-  // offset must be zero for growing axes
+  // TODO this seems to always take the path for growing axes, even if Axes is vector
+  // of variant and types actually held are not growing axes?
+  // index offset must be zero for growing axes
   using IsGrowing = has_growing_axis<Axes>;
   std::fill(indices, indices + size, IsGrowing::value ? 0 : offset);
   for_each_axis(axes, [&, stride = static_cast<std::size_t>(1),

include/boost/histogram/detail/linearize.hpp

@@ -59,7 +59,12 @@ std::size_t linearize(Index& out, const std::size_t stride, const Axis& ax,
   return linearize(opts, out, stride, ax.size(), axis::traits::index(ax, v));
 }
 
-// initial offset of out must be zero
+/**
+  Must be used when axis is potentially growing. Also works for non-growing axis.
+
+  Initial offset of `out` must be zero. We cannot assert on this, because we do not
+  know if this is the first call of `linearize_growth`.
+*/
 template <class Index, class Axis, class Value>
 std::size_t linearize_growth(Index& out, axis::index_type& shift,
                              const std::size_t stride, Axis& a, const Value& v) {

test/CMakeLists.txt

@@ -89,6 +89,7 @@ boost_test(TYPE run SOURCES indexed_test.cpp)
 boost_test(TYPE run SOURCES storage_adaptor_test.cpp)
 boost_test(TYPE run SOURCES unlimited_storage_test.cpp)
 boost_test(TYPE run SOURCES utility_test.cpp)
+boost_test(TYPE run SOURCES issue_327_test.cpp)
 
 find_package(Threads)
 if (Threads_FOUND)

test/Jamfile

@@ -87,6 +87,7 @@ alias cxx14 :
     [ run storage_adaptor_test.cpp ]
     [ run unlimited_storage_test.cpp ]
     [ run utility_test.cpp ]
+    [ run issue_327_test.cpp ]
     ;
 
 alias cxx17 :

test/histogram_fill_test.cpp

@@ -14,6 +14,8 @@
 #include <boost/histogram/axis/category.hpp>
 #include <boost/histogram/axis/integer.hpp>
 #include <boost/histogram/axis/ostream.hpp>
+#include <boost/histogram/axis/regular.hpp>
+#include <boost/histogram/axis/variant.hpp>
 #include <boost/histogram/histogram.hpp>
 #include <boost/histogram/literals.hpp>
 #include <boost/histogram/make_histogram.hpp>
@@ -343,6 +345,61 @@ void run_tests(const std::vector<int>& x, const std::vector<int>& y,
   }
 }
 
+void special_tests() {
+  // 1D growing (bug?)
+  {
+    using axis_type_1 =
+        axis::regular<double, use_default, use_default, axis::option::bitset<11u>>;
+    using axis_type_2 = axis::regular<double>;
+    using axis_type_3 = axis::integer<int>;
+    using axis_type_4 = axis::category<int>;
+    using axis_type_5 = axis::category<std::string>;
+
+    using axis_variant_type =
+        axis::variant<axis_type_1, axis_type_2, axis_type_3, axis_type_4, axis_type_5>;
+
+    auto axes = std::vector<axis_variant_type>({axis_type_1(10, 0, 1)});
+    auto h = histogram<decltype(axes), dense_storage<double>>(axes);
+    auto h2 = h;
+
+    std::vector<int> f1({2});
+    std::vector<int> f2({-1});
+
+    h(2);
+    h(-1);
+
+    h2.fill(f1);
+    h2.fill(f2);
+
+    BOOST_TEST_EQ(h, h2);
+    BOOST_TEST_EQ(sum(h2), 2);
+  }
+
+  // 1D growing (bug?)
+  {
+    using axis_type_1 =
+        axis::regular<double, use_default, use_default, axis::option::bitset<11u>>;
+
+    using axis_variant_type = axis::variant<axis_type_1>;
+
+    auto axes = std::vector<axis_variant_type>({axis_type_1(10, 0, 1)});
+    auto h = histogram<decltype(axes), dense_storage<double>>(axes);
+    auto h2 = h;
+
+    std::vector<int> f1({2});
+    std::vector<int> f2({-1});
+
+    h(2);
+    h(-1);
+
+    h2.fill(f1);
+    h2.fill(f2);
+
+    BOOST_TEST_EQ(h, h2);
+    BOOST_TEST_EQ(sum(h2), 2);
+  }
+}
+
 int main() {
   std::mt19937 gen(1);
   std::normal_distribution<> id(0, 2);
@@ -357,5 +414,7 @@ int main() {
   run_tests<static_tag>(x, y, w);
   run_tests<dynamic_tag>(x, y, w);
 
+  special_tests();
+
   return boost::report_errors();
 }

test/issue_327_test.cpp

@@ -0,0 +1,26 @@
+#include <boost/core/lightweight_test.hpp>
+#include <boost/histogram.hpp>
+#include <boost/histogram/ostream.hpp>
+#include <vector>
+#include "throw_exception.hpp"
+
+namespace bh = boost::histogram;
+using uogrowth_t = decltype(bh::axis::option::growth | bh::axis::option::underflow |
+                            bh::axis::option::overflow);
+
+using arg_t = boost::variant2::variant<std::vector<int>, int>;
+
+int main() {
+  using axis_type =
+      bh::axis::regular<double, bh::use_default, bh::use_default, uogrowth_t>;
+  using axis_variant = bh::axis::variant<axis_type>;
+
+  auto axes = std::vector<axis_variant>({axis_type(10, 0, 1)});
+  auto h = bh::make_histogram_with(std::vector<int>(), axes);
+  BOOST_TEST_EQ(h.rank(), 1);
+
+  std::vector<arg_t> vargs = {-1};
+  h.fill(vargs); // CRASH, using h.fill(-1) or h.fill(args) does not crash.
+
+  return boost::report_errors();
+}