boost/histogram/detail/fill.hpp

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