wip

benchmark/histogram_filling.cpp

@@ -38,17 +38,41 @@ using reg_closed =
 
 class reg_closed_unsafe {
 public:
-  reg_closed_unsafe(axis::index_type n, double start, double stop)
+  using value_type = double;
+
+  reg_closed_unsafe(axis::index_type n, value_type start, value_type stop)
       : min_{start}, delta_{stop - start}, size_{n} {}
 
-  axis::index_type index(double x) const noexcept {
-    // Runs in hot loop, please measure impact of changes
-    auto z = (x - min_) / delta_;
+  axis::index_type index(value_type x) const noexcept {
+    const auto a = 1 / delta_;
+    const auto b = -min_ / delta_;
+    auto z = a * x + b;
     // assume that z < 0 and z > 1 never happens, promised by inclusive()
     if (z == 1) return size() - 1;
     return static_cast<axis::index_type>(z * size());
   }
 
+  template <class T, class I>
+  I* index_transform(T const* begin, T const* end, I* ibegin) const {
+    // constexpr auto n = 128;
+    // value_type buffer[n];
+    // while (begin != end) {
+    //   // write hot computation fused-multiply-add friendly
+    //   const auto a = 1 / delta_;
+    //   const auto b = -min_ / delta_;
+    //   auto vend = buffer;
+    //   // for (; begin != std::min(end, begin + n); ++begin, ++vend) *vend = a * *begin
+    //   +
+    //       // b; for (auto vit = buffer; vit != vend; ++vit, ++ibegin)
+    //       //   *ibegin = std::min(static_cast<axis::index_type>(*vit), size() - 1);
+    //       ++begin;
+    //   *ibegin++ = 0;
+    // }
+    // return ibegin;
+    for (; begin != end; ++begin, ++ibegin) *ibegin = index(*begin);
+    return ibegin;
+  }
+
   axis::index_type size() const noexcept { return size_; }
 
   static constexpr bool inclusive() { return true; }

include/boost/histogram/axis/traits.hpp

@@ -106,11 +106,13 @@ struct variant_access {
   }
 };
 
+// only used if Axis method exists with matching signature
 template <class A>
 decltype(auto) metadata_impl(A&& a, decltype(a.metadata(), 0)) {
   return std::forward<A>(a).metadata();
 }
 
+// fallback option
 template <class A>
 axis::null_type& metadata_impl(A&&, float) {
   static axis::null_type null_value;
@@ -497,6 +499,37 @@ Result width_as(const Axis& axis, index_type index) {
 
 } // namespace traits
 } // namespace axis
+
+namespace detail {
+
+// only used if Axis method exists with matching signature
+template <class A, class VP, class IP>
+IP index_transform_impl(const A& a, VP begin, VP end, IP ibegin,
+                        decltype(a.index_transform(begin, end, ibegin), 0)) {
+  return a.index_transform(begin, end, ibegin);
+}
+
+// fallback option
+template <class A, class VP, class IP>
+IP index_transform_impl(const A& a, VP begin, VP end, IP ibegin, float) {
+  for (auto it = begin; it != end; ++it) *ibegin++ = axis::traits::index(a, *it);
+  return ibegin;
+}
+
+} // namespace detail
+
+namespace axis {
+namespace traits {
+
+template <class Axis, class T, class IndexType>
+IndexType* index_transform(const Axis& axis, const T* begin, const T* end,
+                           IndexType* ibegin) noexcept {
+  return detail::index_transform_impl(axis, begin, end, ibegin, 0);
+}
+
+} // namespace traits
+} // namespace axis
+
 } // namespace histogram
 } // namespace boost
 

include/boost/histogram/detail/fill_n.hpp

@@ -91,15 +91,23 @@ struct index_visitor {
     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
+  template <class T, class IsGrowingOrNot>
+  void call_1(std::false_type, IsGrowingOrNot, const T& iterable) const {
+    // argument is iterable, but value type does not match; fill N values
     const auto* tp = dtl::data(iterable) + start_;
-    for (auto it = begin_; it != begin_ + size_; ++it) call_2(IsGrowing{}, it, *tp++);
+    for (auto it = begin_; it != begin_ + size_; ++it)
+      call_2(IsGrowingOrNot{}, it, *tp++);
   }
 
   template <class T>
-  void call_1(std::true_type, const T& value) const {
+  void call_1(std::false_type, std::false_type, const T& iterable) const {
+    // argument is iterable and axis is not growing; fill N values
+    const auto* tp = dtl::data(iterable) + start_;
+    axis::traits::index_transform(axis_, tp, tp + size_, begin_);
+  }
+
+  template <class T, class IsGrowingOrNot>
+  void call_1(std::true_type, IsGrowingOrNot, 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
@@ -107,7 +115,7 @@ struct index_visitor {
     // 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);
+    call_2(IsGrowingOrNot{}, begin_, value);
     if (is_valid(*begin_)) {
       // since index can be std::size_t or optional_index, must do conversion here
       const auto delta =
@@ -121,7 +129,7 @@ struct index_visitor {
   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);
+           IsGrowing{}, iterable_or_value);
   }
 };
 

test/axis_traits_test.cpp

@@ -30,6 +30,21 @@ struct value_type_deducer<ValueTypeOverride> {
 } // namespace histogram
 } // namespace boost
 
+struct RegularWithoutTransform {
+  index_type index(double x) const { return static_cast<index_type>(x); }
+  index_type size() const { return 10; }
+};
+
+struct RegularWithTransform : RegularWithoutTransform {
+  mutable bool transform_was_called = false;
+  index_type* index_transform(double const* begin, double const* end,
+                              index_type* ibegin) const {
+    transform_was_called = true;
+    for (; begin != end; ++begin, ++ibegin) *ibegin = index(*begin);
+    return ibegin;
+  }
+};
+
 int main() {
   // value_type
   {
@@ -242,6 +257,29 @@ int main() {
     BOOST_TEST_EQ(traits::index(E{0, 3}, "2"), 2);
   }
 
+  // index_transform
+  {
+    double v[3] = {1, 2, 3};
+    index_type i[3];
+
+    const auto iend1 = traits::index_transform(RegularWithoutTransform{}, v, v + 3, i);
+
+    BOOST_TEST_EQ(iend1, i + 3);
+    BOOST_TEST_EQ(i[0], 1);
+    BOOST_TEST_EQ(i[1], 2);
+    BOOST_TEST_EQ(i[2], 3);
+
+    std::fill(i, i + 3, 0);
+
+    RegularWithTransform a;
+    const auto iend2 = traits::index_transform(a, v, v + 3, i);
+    BOOST_TEST_EQ(iend2, i + 3);
+    BOOST_TEST_EQ(i[0], 1);
+    BOOST_TEST_EQ(i[1], 2);
+    BOOST_TEST_EQ(i[2], 3);
+    BOOST_TEST(a.transform_was_called);
+  }
+
   // update
   {
     auto a = integer<int, null_type, option::growth_t>();