histogram/detail/fill.hpp

0001 // Copyright 2015-2018 Hans Dembinski
0002 //
0003 // Distributed under the Boost Software License, Version 1.0.
0004 // (See accompanying file LICENSE_1_0.txt
0005 // or copy at http://www.boost.org/LICENSE_1_0.txt)
0006
0007 #ifndef BOOST_HISTOGRAM_DETAIL_FILL_HPP
0008 #define BOOST_HISTOGRAM_DETAIL_FILL_HPP
0009
0010 #include <algorithm>
0011 #include <boost/config/workaround.hpp>
0012 #include <boost/histogram/axis/traits.hpp>
0013 #include <boost/histogram/axis/variant.hpp>
0014 #include <boost/histogram/detail/argument_traits.hpp>
0015 #include <boost/histogram/detail/axes.hpp>
0016 #include <boost/histogram/detail/linearize.hpp>
0017 #include <boost/histogram/detail/make_default.hpp>
0018 #include <boost/histogram/detail/optional_index.hpp>
0019 #include <boost/histogram/detail/priority.hpp>
0020 #include <boost/histogram/detail/tuple_slice.hpp>
0021 #include <boost/histogram/fwd.hpp>
0022 #include <boost/mp11/algorithm.hpp>
0023 #include <boost/mp11/integral.hpp>
0024 #include <boost/mp11/tuple.hpp>
0025 #include <boost/mp11/utility.hpp>
0026 #include <cassert>
0027 #include <mutex>
0028 #include <tuple>
0029 #include <type_traits>
0030
0031 namespace boost {
0032 namespace histogram {
0033 namespace detail {
0034
0035 template <class T, class U>
0036 struct sample_args_passed_vs_expected;
0037
0038 template <class... Passed, class... Expected>
0039 struct sample_args_passed_vs_expected<std::tuple<Passed...>, std::tuple<Expected...>> {
0040   static_assert(!(sizeof...(Expected) > 0 && sizeof...(Passed) == 0),
0041                 "error: accumulator requires samples, but sample argument is missing");
0042   static_assert(
0043       !(sizeof...(Passed) > 0 && sizeof...(Expected) == 0),
0044       "error: accumulator does not accept samples, but sample argument is passed");
0045   static_assert(sizeof...(Passed) == sizeof...(Expected),
0046                 "error: numbers of passed and expected sample arguments differ");
0047   static_assert(
0048       std::is_convertible<std::tuple<Passed...>, std::tuple<Expected...>>::value,
0049       "error: sample argument(s) not convertible to accumulator argument(s)");
0050 };
0051
0052 template <class A>
0053 struct storage_grower {
0054   const A& axes_;
0055   struct {
0056     axis::index_type idx, old_extent;
0057     std::size_t new_stride;
0058   } data_[buffer_size<A>::value];
0059   std::size_t new_size_;
0060
0061   storage_grower(const A& axes) noexcept : axes_(axes) {}
0062
0063   void from_shifts(const axis::index_type* shifts) noexcept {
0064     auto dit = data_;
0065     std::size_t s = 1;
0066     for_each_axis(axes_, [&](const auto& a) {
0067       const auto n = axis::traits::extent(a);
0068       *dit++ = {0, n - std::abs(*shifts++), s};
0069       s *= n;
0070     });
0071     new_size_ = s;
0072   }
0073
0074   // must be extents before any shifts were applied
0075   void from_extents(const axis::index_type* old_extents) noexcept {
0076     auto dit = data_;
0077     std::size_t s = 1;
0078     for_each_axis(axes_, [&](const auto& a) {
0079       const auto n = axis::traits::extent(a);
0080       *dit++ = {0, *old_extents++, s};
0081       s *= n;
0082     });
0083     new_size_ = s;
0084   }
0085
0086   template <class S>
0087   void apply(S& storage, const axis::index_type* shifts) {
0088     auto new_storage = make_default(storage);
0089     new_storage.reset(new_size_);
0090     const auto dlast = data_ + axes_rank(axes_) - 1;
0091     for (auto&& x : storage) {
0092       auto ns = new_storage.begin();
0093       auto sit = shifts;
0094       auto dit = data_;
0095       for_each_axis(axes_, [&](const auto& a) {
0096         using opt = axis::traits::get_options<std::decay_t<decltype(a)>>;
0097         if (opt::test(axis::option::underflow)) {
0098           if (dit->idx == 0) {
0099             // axis has underflow and we are in the underflow bin:
0100             // keep storage pointer unchanged
0101             ++dit;
0102             ++sit;
0103             return;
0104           }
0105         }
0106         if (opt::test(axis::option::overflow)) {
0107           if (dit->idx == dit->old_extent - 1) {
0108             // axis has overflow and we are in the overflow bin:
0109             // move storage pointer to corresponding overflow bin position
0110             ns += (axis::traits::extent(a) - 1) * dit->new_stride;
0111             ++dit;
0112             ++sit;
0113             return;
0114           }
0115         }
0116         // we are in a normal bin:
0117         // move storage pointer to index position; apply positive shifts if any
0118         ns += (dit->idx + (*sit >= 0 ? *sit : 0)) * dit->new_stride;
0119         ++dit;
0120         ++sit;
0121       });
0122       // assign old value to new location
0123       *ns = x;
0124       // advance multi-dimensional index
0125       dit = data_;
0126       ++dit->idx;
0127       while (dit != dlast && dit->idx == dit->old_extent) {
0128         dit->idx = 0;
0129         ++(++dit)->idx;
0130       }
0131     }
0132     storage = std::move(new_storage);
0133   }
0134 };
0135
0136 template <class T, class... Us>
0137 auto fill_storage_element_impl(priority<2>, T&& t, const Us&... args) noexcept
0138     -> decltype(t(args...), void()) {
0139   t(args...);
0140 }
0141
0142 template <class T, class U>
0143 auto fill_storage_element_impl(priority<1>, T&& t, const weight_type<U>& w) noexcept
0144     -> decltype(t += w, void()) {
0145   t += w;
0146 }
0147
0148 // fallback for arithmetic types and accumulators that do not handle the weight
0149 template <class T, class U>
0150 auto fill_storage_element_impl(priority<0>, T&& t, const weight_type<U>& w) noexcept
0151     -> decltype(t += w.value, void()) {
0152   t += w.value;
0153 }
0154
0155 template <class T>
0156 auto fill_storage_element_impl(priority<1>, T&& t) noexcept -> decltype(++t, void()) {
0157   ++t;
0158 }
0159
0160 template <class T, class... Us>
0161 void fill_storage_element(T&& t, const Us&... args) noexcept {
0162   fill_storage_element_impl(priority<2>{}, std::forward<T>(t), args...);
0163 }
0164
0165 // t may be a proxy and then it is an rvalue reference, not an lvalue reference
0166 template <class IW, class IS, class T, class U>
0167 void fill_storage_2(IW, IS, T&& t, U&& u) noexcept {
0168   mp11::tuple_apply(
0169       [&](const auto&... args) {
0170         fill_storage_element(std::forward<T>(t), std::get<IW::value>(u), args...);
0171       },
0172       std::get<IS::value>(u).value);
0173 }
0174
0175 // t may be a proxy and then it is an rvalue reference, not an lvalue reference
0176 template <class IS, class T, class U>
0177 void fill_storage_2(mp11::mp_int<-1>, IS, T&& t, const U& u) noexcept {
0178   mp11::tuple_apply(
0179       [&](const auto&... args) { fill_storage_element(std::forward<T>(t), args...); },
0180       std::get<IS::value>(u).value);
0181 }
0182
0183 // t may be a proxy and then it is an rvalue reference, not an lvalue reference
0184 template <class IW, class T, class U>
0185 void fill_storage_2(IW, mp11::mp_int<-1>, T&& t, const U& u) noexcept {
0186   fill_storage_element(std::forward<T>(t), std::get<IW::value>(u));
0187 }
0188
0189 // t may be a proxy and then it is an rvalue reference, not an lvalue reference
0190 template <class T, class U>
0191 void fill_storage_2(mp11::mp_int<-1>, mp11::mp_int<-1>, T&& t, const U&) noexcept {
0192   fill_storage_element(std::forward<T>(t));
0193 }
0194
0195 template <class IW, class IS, class Storage, class Index, class Args>
0196 auto fill_storage(IW, IS, Storage& s, const Index idx, const Args& a) noexcept {
0197   if (is_valid(idx)) {
0198     assert(idx < s.size());
0199     fill_storage_2(IW{}, IS{}, s[idx], a);
0200     return s.begin() + idx;
0201   }
0202   return s.end();
0203 }
0204
0205 template <int S, int N>
0206 struct linearize_args {
0207   template <class Index, class A, class Args>
0208   static void impl(mp11::mp_int<N>, Index&, const std::size_t, A&, const Args&) {}
0209
0210   template <int I, class Index, class A, class Args>
0211   static void impl(mp11::mp_int<I>, Index& o, const std::size_t s, A& ax,
0212                    const Args& args) {
0213     const auto e = linearize(o, s, axis_get<I>(ax), std::get<(S + I)>(args));
0214     impl(mp11::mp_int<(I + 1)>{}, o, s * e, ax, args);
0215   }
0216
0217   template <class Index, class A, class Args>
0218   static void apply(Index& o, A& ax, const Args& args) {
0219     impl(mp11::mp_int<0>{}, o, 1, ax, args);
0220   }
0221 };
0222
0223 template <int S>
0224 struct linearize_args<S, 1> {
0225   template <class Index, class A, class Args>
0226   static void apply(Index& o, A& ax, const Args& args) {
0227     linearize(o, 1, axis_get<0>(ax), std::get<S>(args));
0228   }
0229 };
0230
0231 template <class A>
0232 constexpr unsigned(min)(const unsigned n) noexcept {
0233   constexpr unsigned a = buffer_size<A>::value;
0234   return a < n ? a : n;
0235 }
0236
0237 // not growing
0238 template <class ArgTraits, class Storage, class Axes, class Args>
0239 auto fill_2(ArgTraits, mp11::mp_false, const std::size_t offset, Storage& st,
0240             const Axes& axes, const Args& args) {
0241   mp11::mp_if<has_non_inclusive_axis<Axes>, optional_index, std::size_t> idx{offset};
0242   linearize_args<ArgTraits::start::value, min<Axes>(ArgTraits::nargs::value)>::apply(
0243       idx, axes, args);
0244   return fill_storage(typename ArgTraits::wpos{}, typename ArgTraits::spos{}, st, idx,
0245                       args);
0246 }
0247
0248 // at least one axis is growing
0249 template <class ArgTraits, class Storage, class Axes, class Args>
0250 auto fill_2(ArgTraits, mp11::mp_true, const std::size_t, Storage& st, Axes& axes,
0251             const Args& args) {
0252   std::array<axis::index_type, ArgTraits::nargs::value> shifts;
0253   // offset must be zero for linearize_growth (value of offset argument is ignored)
0254   mp11::mp_if<has_non_inclusive_axis<Axes>, optional_index, std::size_t> idx{0};
0255   std::size_t stride = 1;
0256   bool update_needed = false;
0257   mp11::mp_for_each<mp11::mp_iota_c<min<Axes>(ArgTraits::nargs::value)>>([&](auto i) {
0258     auto& ax = axis_get<i>(axes);
0259     const auto extent = linearize_growth(idx, shifts[i], stride, ax,
0260                                          std::get<(ArgTraits::start::value + i)>(args));
0261     update_needed |= shifts[i] != 0;
0262     stride *= extent;
0263   });
0264   if (update_needed) {
0265     storage_grower<Axes> g(axes);
0266     g.from_shifts(shifts.data());
0267     g.apply(st, shifts.data());
0268   }
0269   return fill_storage(typename ArgTraits::wpos{}, typename ArgTraits::spos{}, st, idx,
0270                       args);
0271 }
0272
0273 // pack original args tuple into another tuple (which is unpacked later)
0274 template <int Start, int Size, class IW, class IS, class Args>
0275 decltype(auto) pack_args(IW, IS, const Args& args) noexcept {
0276   return std::make_tuple(tuple_slice<Start, Size>(args), std::get<IW::value>(args),
0277                          std::get<IS::value>(args));
0278 }
0279
0280 template <int Start, int Size, class IW, class Args>
0281 decltype(auto) pack_args(IW, mp11::mp_int<-1>, const Args& args) noexcept {
0282   return std::make_tuple(tuple_slice<Start, Size>(args), std::get<IW::value>(args));
0283 }
0284
0285 template <int Start, int Size, class IS, class Args>
0286 decltype(auto) pack_args(mp11::mp_int<-1>, IS, const Args& args) noexcept {
0287   return std::make_tuple(tuple_slice<Start, Size>(args), std::get<IS::value>(args));
0288 }
0289
0290 template <int Start, int Size, class Args>
0291 decltype(auto) pack_args(mp11::mp_int<-1>, mp11::mp_int<-1>, const Args& args) noexcept {
0292   return std::make_tuple(args);
0293 }
0294
0295 #if BOOST_WORKAROUND(BOOST_MSVC, >= 0)
0296 #pragma warning(disable : 4702) // fixing warning would reduce code readability a lot
0297 #endif
0298
0299 template <class ArgTraits, class S, class A, class Args>
0300 auto fill(std::true_type, ArgTraits, const std::size_t offset, S& storage, A& axes,
0301           const Args& args) -> typename S::iterator {
0302   using growing = has_growing_axis<A>;
0303
0304   // Sometimes we need to pack the tuple into another tuple:
0305   // - histogram contains one axis which accepts tuple
0306   // - user passes tuple to fill(...)
0307   // Tuple is normally unpacked and arguments are processed, this causes pos::nargs > 1.
0308   // Now we pack tuple into another tuple so that original tuple is send to axis.
0309   // Notes:
0310   // - has nice side-effect of making histogram::operator(1, 2) work as well
0311   // - cannot detect call signature of axis at compile-time in all configurations
0312   //   (axis::variant provides generic call interface and hides concrete
0313   //   interface), so we throw at runtime if incompatible argument is passed (e.g.
0314   //   3d tuple)
0315
0316   if (axes_rank(axes) == ArgTraits::nargs::value)
0317     return fill_2(ArgTraits{}, growing{}, offset, storage, axes, args);
0318   else if (axes_rank(axes) == 1 &&
0319            axis::traits::rank(axis_get<0>(axes)) == ArgTraits::nargs::value)
0320     return fill_2(
0321         argument_traits_holder<
0322             1, 0, (ArgTraits::wpos::value >= 0 ? 1 : -1),
0323             (ArgTraits::spos::value >= 0 ? (ArgTraits::wpos::value >= 0 ? 2 : 1) : -1),
0324             typename ArgTraits::sargs>{},
0325         growing{}, offset, storage, axes,
0326         pack_args<ArgTraits::start::value, ArgTraits::nargs::value>(
0327             typename ArgTraits::wpos{}, typename ArgTraits::spos{}, args));
0328   return BOOST_THROW_EXCEPTION(
0329              std::invalid_argument("number of arguments != histogram rank")),
0330          storage.end();
0331 }
0332
0333 #if BOOST_WORKAROUND(BOOST_MSVC, >= 0)
0334 #pragma warning(default : 4702)
0335 #endif
0336
0337 // empty implementation for bad arguments to stop compiler from showing internals
0338 template <class ArgTraits, class S, class A, class Args>
0339 auto fill(std::false_type, ArgTraits, const std::size_t, S& storage, A&, const Args&) ->
0340     typename S::iterator {
0341   return storage.end();
0342 }
0343
0344 } // namespace detail
0345 } // namespace histogram
0346 } // namespace boost
0347
0348 #endif