fill.hpp 12 KB

123456789101112131415161718192021222324252627282930313233343536373839404142434445464748495051525354555657585960616263646566676869707172737475767778798081828384858687888990919293949596979899100101102103104105106107108109110111112113114115116117118119120121122123124125126127128129130131132133134135136137138139140141142143144145146147148149150151152153154155156157158159160161162163164165166167168169170171172173174175176177178179180181182183184185186187188189190191192193194195196197198199200201202203204205206207208209210211212213214215216217218219220221222223224225226227228229230231232233234235236237238239240241242243244245246247248249250251252253254255256257258259260261262263264265266267268269270271272273274275276277278279280281282283284285286287288289290291292293294295296297298299300301302303304305306307308309310311312313314315316317318319320321322323324325326327328329330331332333334335336337338339340341342343344345346347348
  1. // Copyright 2015-2018 Hans Dembinski
  2. //
  3. // Distributed under the Boost Software License, Version 1.0.
  4. // (See accompanying file LICENSE_1_0.txt
  5. // or copy at http://www.boost.org/LICENSE_1_0.txt)
  6. #ifndef BOOST_HISTOGRAM_DETAIL_FILL_HPP
  7. #define BOOST_HISTOGRAM_DETAIL_FILL_HPP
  8. #include <algorithm>
  9. #include <boost/config/workaround.hpp>
  10. #include <boost/histogram/axis/traits.hpp>
  11. #include <boost/histogram/axis/variant.hpp>
  12. #include <boost/histogram/detail/argument_traits.hpp>
  13. #include <boost/histogram/detail/axes.hpp>
  14. #include <boost/histogram/detail/linearize.hpp>
  15. #include <boost/histogram/detail/make_default.hpp>
  16. #include <boost/histogram/detail/optional_index.hpp>
  17. #include <boost/histogram/detail/priority.hpp>
  18. #include <boost/histogram/detail/tuple_slice.hpp>
  19. #include <boost/histogram/fwd.hpp>
  20. #include <boost/mp11/algorithm.hpp>
  21. #include <boost/mp11/integral.hpp>
  22. #include <boost/mp11/tuple.hpp>
  23. #include <boost/mp11/utility.hpp>
  24. #include <cassert>
  25. #include <mutex>
  26. #include <tuple>
  27. #include <type_traits>
  28. namespace boost {
  29. namespace histogram {
  30. namespace detail {
  31. template <class T, class U>
  32. struct sample_args_passed_vs_expected;
  33. template <class... Passed, class... Expected>
  34. struct sample_args_passed_vs_expected<std::tuple<Passed...>, std::tuple<Expected...>> {
  35. static_assert(!(sizeof...(Expected) > 0 && sizeof...(Passed) == 0),
  36. "error: accumulator requires samples, but sample argument is missing");
  37. static_assert(
  38. !(sizeof...(Passed) > 0 && sizeof...(Expected) == 0),
  39. "error: accumulator does not accept samples, but sample argument is passed");
  40. static_assert(sizeof...(Passed) == sizeof...(Expected),
  41. "error: numbers of passed and expected sample arguments differ");
  42. static_assert(
  43. std::is_convertible<std::tuple<Passed...>, std::tuple<Expected...>>::value,
  44. "error: sample argument(s) not convertible to accumulator argument(s)");
  45. };
  46. template <class A>
  47. struct storage_grower {
  48. const A& axes_;
  49. struct {
  50. axis::index_type idx, old_extent;
  51. std::size_t new_stride;
  52. } data_[buffer_size<A>::value];
  53. std::size_t new_size_;
  54. storage_grower(const A& axes) noexcept : axes_(axes) {}
  55. void from_shifts(const axis::index_type* shifts) noexcept {
  56. auto dit = data_;
  57. std::size_t s = 1;
  58. for_each_axis(axes_, [&](const auto& a) {
  59. const auto n = axis::traits::extent(a);
  60. *dit++ = {0, n - std::abs(*shifts++), s};
  61. s *= n;
  62. });
  63. new_size_ = s;
  64. }
  65. // must be extents before any shifts were applied
  66. void from_extents(const axis::index_type* old_extents) noexcept {
  67. auto dit = data_;
  68. std::size_t s = 1;
  69. for_each_axis(axes_, [&](const auto& a) {
  70. const auto n = axis::traits::extent(a);
  71. *dit++ = {0, *old_extents++, s};
  72. s *= n;
  73. });
  74. new_size_ = s;
  75. }
  76. template <class S>
  77. void apply(S& storage, const axis::index_type* shifts) {
  78. auto new_storage = make_default(storage);
  79. new_storage.reset(new_size_);
  80. const auto dlast = data_ + axes_rank(axes_) - 1;
  81. for (auto&& x : storage) {
  82. auto ns = new_storage.begin();
  83. auto sit = shifts;
  84. auto dit = data_;
  85. for_each_axis(axes_, [&](const auto& a) {
  86. using opt = axis::traits::get_options<std::decay_t<decltype(a)>>;
  87. if (opt::test(axis::option::underflow)) {
  88. if (dit->idx == 0) {
  89. // axis has underflow and we are in the underflow bin:
  90. // keep storage pointer unchanged
  91. ++dit;
  92. ++sit;
  93. return;
  94. }
  95. }
  96. if (opt::test(axis::option::overflow)) {
  97. if (dit->idx == dit->old_extent - 1) {
  98. // axis has overflow and we are in the overflow bin:
  99. // move storage pointer to corresponding overflow bin position
  100. ns += (axis::traits::extent(a) - 1) * dit->new_stride;
  101. ++dit;
  102. ++sit;
  103. return;
  104. }
  105. }
  106. // we are in a normal bin:
  107. // move storage pointer to index position; apply positive shifts if any
  108. ns += (dit->idx + (*sit >= 0 ? *sit : 0)) * dit->new_stride;
  109. ++dit;
  110. ++sit;
  111. });
  112. // assign old value to new location
  113. *ns = x;
  114. // advance multi-dimensional index
  115. dit = data_;
  116. ++dit->idx;
  117. while (dit != dlast && dit->idx == dit->old_extent) {
  118. dit->idx = 0;
  119. ++(++dit)->idx;
  120. }
  121. }
  122. storage = std::move(new_storage);
  123. }
  124. };
  125. template <class T, class... Us>
  126. auto fill_storage_element_impl(priority<2>, T&& t, const Us&... args) noexcept
  127. -> decltype(t(args...), void()) {
  128. t(args...);
  129. }
  130. template <class T, class U>
  131. auto fill_storage_element_impl(priority<1>, T&& t, const weight_type<U>& w) noexcept
  132. -> decltype(t += w, void()) {
  133. t += w;
  134. }
  135. // fallback for arithmetic types and accumulators that do not handle the weight
  136. template <class T, class U>
  137. auto fill_storage_element_impl(priority<0>, T&& t, const weight_type<U>& w) noexcept
  138. -> decltype(t += w.value, void()) {
  139. t += w.value;
  140. }
  141. template <class T>
  142. auto fill_storage_element_impl(priority<1>, T&& t) noexcept -> decltype(++t, void()) {
  143. ++t;
  144. }
  145. template <class T, class... Us>
  146. void fill_storage_element(T&& t, const Us&... args) noexcept {
  147. fill_storage_element_impl(priority<2>{}, std::forward<T>(t), args...);
  148. }
  149. // t may be a proxy and then it is an rvalue reference, not an lvalue reference
  150. template <class IW, class IS, class T, class U>
  151. void fill_storage_2(IW, IS, T&& t, U&& u) noexcept {
  152. mp11::tuple_apply(
  153. [&](const auto&... args) {
  154. fill_storage_element(std::forward<T>(t), std::get<IW::value>(u), args...);
  155. },
  156. std::get<IS::value>(u).value);
  157. }
  158. // t may be a proxy and then it is an rvalue reference, not an lvalue reference
  159. template <class IS, class T, class U>
  160. void fill_storage_2(mp11::mp_int<-1>, IS, T&& t, const U& u) noexcept {
  161. mp11::tuple_apply(
  162. [&](const auto&... args) { fill_storage_element(std::forward<T>(t), args...); },
  163. std::get<IS::value>(u).value);
  164. }
  165. // t may be a proxy and then it is an rvalue reference, not an lvalue reference
  166. template <class IW, class T, class U>
  167. void fill_storage_2(IW, mp11::mp_int<-1>, T&& t, const U& u) noexcept {
  168. fill_storage_element(std::forward<T>(t), std::get<IW::value>(u));
  169. }
  170. // t may be a proxy and then it is an rvalue reference, not an lvalue reference
  171. template <class T, class U>
  172. void fill_storage_2(mp11::mp_int<-1>, mp11::mp_int<-1>, T&& t, const U&) noexcept {
  173. fill_storage_element(std::forward<T>(t));
  174. }
  175. template <class IW, class IS, class Storage, class Index, class Args>
  176. auto fill_storage(IW, IS, Storage& s, const Index idx, const Args& a) noexcept {
  177. if (is_valid(idx)) {
  178. assert(idx < s.size());
  179. fill_storage_2(IW{}, IS{}, s[idx], a);
  180. return s.begin() + idx;
  181. }
  182. return s.end();
  183. }
  184. template <int S, int N>
  185. struct linearize_args {
  186. template <class Index, class A, class Args>
  187. static void impl(mp11::mp_int<N>, Index&, const std::size_t, A&, const Args&) {}
  188. template <int I, class Index, class A, class Args>
  189. static void impl(mp11::mp_int<I>, Index& o, const std::size_t s, A& ax,
  190. const Args& args) {
  191. const auto e = linearize(o, s, axis_get<I>(ax), std::get<(S + I)>(args));
  192. impl(mp11::mp_int<(I + 1)>{}, o, s * e, ax, args);
  193. }
  194. template <class Index, class A, class Args>
  195. static void apply(Index& o, A& ax, const Args& args) {
  196. impl(mp11::mp_int<0>{}, o, 1, ax, args);
  197. }
  198. };
  199. template <int S>
  200. struct linearize_args<S, 1> {
  201. template <class Index, class A, class Args>
  202. static void apply(Index& o, A& ax, const Args& args) {
  203. linearize(o, 1, axis_get<0>(ax), std::get<S>(args));
  204. }
  205. };
  206. template <class A>
  207. constexpr unsigned min(const unsigned n) noexcept {
  208. constexpr unsigned a = buffer_size<A>::value;
  209. return a < n ? a : n;
  210. }
  211. // not growing
  212. template <class ArgTraits, class Storage, class Axes, class Args>
  213. auto fill_2(ArgTraits, mp11::mp_false, const std::size_t offset, Storage& st,
  214. const Axes& axes, const Args& args) {
  215. mp11::mp_if<has_non_inclusive_axis<Axes>, optional_index, std::size_t> idx{offset};
  216. linearize_args<ArgTraits::start::value, min<Axes>(ArgTraits::nargs::value)>::apply(
  217. idx, axes, args);
  218. return fill_storage(typename ArgTraits::wpos{}, typename ArgTraits::spos{}, st, idx,
  219. args);
  220. }
  221. // at least one axis is growing
  222. template <class ArgTraits, class Storage, class Axes, class Args>
  223. auto fill_2(ArgTraits, mp11::mp_true, const std::size_t, Storage& st, Axes& axes,
  224. const Args& args) {
  225. std::array<axis::index_type, ArgTraits::nargs::value> shifts;
  226. // offset must be zero for linearize_growth
  227. mp11::mp_if<has_non_inclusive_axis<Axes>, optional_index, std::size_t> idx{0};
  228. std::size_t stride = 1;
  229. bool update_needed = false;
  230. mp11::mp_for_each<mp11::mp_iota_c<min<Axes>(ArgTraits::nargs::value)>>([&](auto i) {
  231. auto& ax = axis_get<i>(axes);
  232. const auto extent = linearize_growth(idx, shifts[i], stride, ax,
  233. std::get<(ArgTraits::start::value + i)>(args));
  234. update_needed |= shifts[i] != 0;
  235. stride *= extent;
  236. });
  237. if (update_needed) {
  238. storage_grower<Axes> g(axes);
  239. g.from_shifts(shifts.data());
  240. g.apply(st, shifts.data());
  241. }
  242. return fill_storage(typename ArgTraits::wpos{}, typename ArgTraits::spos{}, st, idx,
  243. args);
  244. }
  245. // pack original args tuple into another tuple (which is unpacked later)
  246. template <int Start, int Size, class IW, class IS, class Args>
  247. decltype(auto) pack_args(IW, IS, const Args& args) noexcept {
  248. return std::make_tuple(tuple_slice<Start, Size>(args), std::get<IW::value>(args),
  249. std::get<IS::value>(args));
  250. }
  251. template <int Start, int Size, class IW, class Args>
  252. decltype(auto) pack_args(IW, mp11::mp_int<-1>, const Args& args) noexcept {
  253. return std::make_tuple(tuple_slice<Start, Size>(args), std::get<IW::value>(args));
  254. }
  255. template <int Start, int Size, class IS, class Args>
  256. decltype(auto) pack_args(mp11::mp_int<-1>, IS, const Args& args) noexcept {
  257. return std::make_tuple(tuple_slice<Start, Size>(args), std::get<IS::value>(args));
  258. }
  259. template <int Start, int Size, class Args>
  260. decltype(auto) pack_args(mp11::mp_int<-1>, mp11::mp_int<-1>, const Args& args) noexcept {
  261. return std::make_tuple(args);
  262. }
  263. #if BOOST_WORKAROUND(BOOST_MSVC, >= 0)
  264. #pragma warning(disable : 4702) // fixing warning would reduce code readability a lot
  265. #endif
  266. template <class ArgTraits, class S, class A, class Args>
  267. auto fill(std::true_type, ArgTraits, const std::size_t offset, S& storage, A& axes,
  268. const Args& args) -> typename S::iterator {
  269. using growing = has_growing_axis<A>;
  270. // Sometimes we need to pack the tuple into another tuple:
  271. // - histogram contains one axis which accepts tuple
  272. // - user passes tuple to fill(...)
  273. // Tuple is normally unpacked and arguments are processed, this causes pos::nargs > 1.
  274. // Now we pack tuple into another tuple so that original tuple is send to axis.
  275. // Notes:
  276. // - has nice side-effect of making histogram::operator(1, 2) work as well
  277. // - cannot detect call signature of axis at compile-time in all configurations
  278. // (axis::variant provides generic call interface and hides concrete
  279. // interface), so we throw at runtime if incompatible argument is passed (e.g.
  280. // 3d tuple)
  281. if (axes_rank(axes) == ArgTraits::nargs::value)
  282. return fill_2(ArgTraits{}, growing{}, offset, storage, axes, args);
  283. else if (axes_rank(axes) == 1 &&
  284. axis::traits::rank(axis_get<0>(axes)) == ArgTraits::nargs::value)
  285. return fill_2(
  286. argument_traits_holder<
  287. 1, 0, (ArgTraits::wpos::value >= 0 ? 1 : -1),
  288. (ArgTraits::spos::value >= 0 ? (ArgTraits::wpos::value >= 0 ? 2 : 1) : -1),
  289. typename ArgTraits::sargs>{},
  290. growing{}, offset, storage, axes,
  291. pack_args<ArgTraits::start::value, ArgTraits::nargs::value>(
  292. typename ArgTraits::wpos{}, typename ArgTraits::spos{}, args));
  293. return BOOST_THROW_EXCEPTION(
  294. std::invalid_argument("number of arguments != histogram rank")),
  295. storage.end();
  296. }
  297. #if BOOST_WORKAROUND(BOOST_MSVC, >= 0)
  298. #pragma warning(default : 4702)
  299. #endif
  300. // empty implementation for bad arguments to stop compiler from showing internals
  301. template <class ArgTraits, class S, class A, class Args>
  302. auto fill(std::false_type, ArgTraits, const std::size_t, S& storage, A&, const Args&) ->
  303. typename S::iterator {
  304. return storage.end();
  305. }
  306. } // namespace detail
  307. } // namespace histogram
  308. } // namespace boost
  309. #endif