boost/histogram/detail/axes.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_AXES_HPP
   8 #define BOOST_HISTOGRAM_DETAIL_AXES_HPP
   9
  10 #include <array>
  11 #include <boost/assert.hpp>
  12 #include <boost/core/nvp.hpp>
  13 #include <boost/histogram/axis/traits.hpp>
  14 #include <boost/histogram/axis/variant.hpp>
  15 #include <boost/histogram/detail/make_default.hpp>
  16 #include <boost/histogram/detail/optional_index.hpp>
  17 #include <boost/histogram/detail/static_if.hpp>
  18 #include <boost/histogram/fwd.hpp>
  19 #include <boost/mp11/algorithm.hpp>
  20 #include <boost/mp11/list.hpp>
  21 #include <boost/mp11/tuple.hpp>
  22 #include <boost/mp11/utility.hpp>
  23 #include <boost/throw_exception.hpp>
  24 #include <stdexcept>
  25 #include <string>
  26 #include <tuple>
  27 #include <type_traits>
  28
  29 /* Most of the histogram code is generic and works for any number of axes. Buffers with a
  30  * fixed maximum capacity are used in some places, which have a size equal to the rank of
  31  * a histogram. The buffers are statically allocated to improve performance, which means
  32  * that they need a preset maximum capacity. 32 seems like a safe upper limit for the rank
  33  * (you can nevertheless increase it here if necessary): the simplest non-trivial axis has
  34  * 2 bins; even if counters are used which need only a byte of storage per bin, this still
  35  * corresponds to 4 GB of storage.
  36  */
  37 #ifndef BOOST_HISTOGRAM_DETAIL_AXES_LIMIT
  38 #define BOOST_HISTOGRAM_DETAIL_AXES_LIMIT 32
  39 #endif
  40
  41 namespace boost {
  42 namespace histogram {
  43 namespace detail {
  44
  45 template <class T>
  46 unsigned axes_rank(const T& axes) {
  47   using std::begin;
  48   using std::end;
  49   return static_cast<unsigned>(std::distance(begin(axes), end(axes)));
  50 }
  51
  52 template <class... Ts>
  53 constexpr unsigned axes_rank(const std::tuple<Ts...>&) {
  54   return static_cast<unsigned>(sizeof...(Ts));
  55 }
  56
  57 template <class T>
  58 void throw_if_axes_is_too_large(const T& axes) {
  59   if (axes_rank(axes) > BOOST_HISTOGRAM_DETAIL_AXES_LIMIT)
  60     BOOST_THROW_EXCEPTION(
  61         std::invalid_argument("length of axis vector exceeds internal buffers, "
  62                               "recompile with "
  63                               "-DBOOST_HISTOGRAM_DETAIL_AXES_LIMIT=<new max size> "
  64                               "to increase internal buffers"));
  65 }
  66
  67 // tuple is never too large because internal buffers adapt to size of tuple
  68 template <class... Ts>
  69 void throw_if_axes_is_too_large(const std::tuple<Ts...>&) {}
  70
  71 template <unsigned N, class... Ts>
  72 decltype(auto) axis_get(std::tuple<Ts...>& axes) {
  73   return std::get<N>(axes);
  74 }
  75
  76 template <unsigned N, class... Ts>
  77 decltype(auto) axis_get(const std::tuple<Ts...>& axes) {
  78   return std::get<N>(axes);
  79 }
  80
  81 template <unsigned N, class T>
  82 decltype(auto) axis_get(T& axes) {
  83   return axes[N];
  84 }
  85
  86 template <unsigned N, class T>
  87 decltype(auto) axis_get(const T& axes) {
  88   return axes[N];
  89 }
  90
  91 template <class... Ts>
  92 auto axis_get(std::tuple<Ts...>& axes, const unsigned i) {
  93   constexpr auto S = sizeof...(Ts);
  94   using V = mp11::mp_unique<axis::variant<Ts*...>>;
  95   return mp11::mp_with_index<S>(i, [&axes](auto i) { return V(&std::get<i>(axes)); });
  96 }
  97
  98 template <class... Ts>
  99 auto axis_get(const std::tuple<Ts...>& axes, const unsigned i) {
 100   constexpr auto S = sizeof...(Ts);
 101   using V = mp11::mp_unique<axis::variant<const Ts*...>>;
 102   return mp11::mp_with_index<S>(i, [&axes](auto i) { return V(&std::get<i>(axes)); });
 103 }
 104
 105 template <class T>
 106 decltype(auto) axis_get(T& axes, const unsigned i) {
 107   return axes[i];
 108 }
 109
 110 template <class T>
 111 decltype(auto) axis_get(const T& axes, const unsigned i) {
 112   return axes[i];
 113 }
 114
 115 template <class... Ts, class... Us>
 116 bool axes_equal(const std::tuple<Ts...>& ts, const std::tuple<Us...>& us) {
 117   using namespace ::boost::mp11;
 118   return static_if<std::is_same<mp_list<Ts...>, mp_list<Us...>>>(
 119       [](const auto& ts, const auto& us) {
 120         using N = mp_size<std::decay_t<decltype(ts)>>;
 121         bool equal = true;
 122         mp_for_each<mp_iota<N>>(
 123             [&](auto I) { equal &= relaxed_equal(std::get<I>(ts), std::get<I>(us)); });
 124         return equal;
 125       },
 126       [](const auto&, const auto&) { return false; }, ts, us);
 127 }
 128
 129 template <class T, class... Us>
 130 bool axes_equal(const T& t, const std::tuple<Us...>& u) {
 131   using namespace ::boost::mp11;
 132   if (t.size() != sizeof...(Us)) return false;
 133   bool equal = true;
 134   mp_for_each<mp_iota_c<sizeof...(Us)>>([&](auto I) { equal &= t[I] == std::get<I>(u); });
 135   return equal;
 136 }
 137
 138 template <class... Ts, class U>
 139 bool axes_equal(const std::tuple<Ts...>& t, const U& u) {
 140   return axes_equal(u, t);
 141 }
 142
 143 template <class T, class U>
 144 bool axes_equal(const T& t, const U& u) {
 145   if (t.size() != u.size()) return false;
 146   return std::equal(t.begin(), t.end(), u.begin());
 147 }
 148
 149 template <class... Ts, class... Us>
 150 void axes_assign(std::tuple<Ts...>& t, const std::tuple<Us...>& u) {
 151   using namespace ::boost::mp11;
 152   static_if<std::is_same<mp_list<Ts...>, mp_list<Us...>>>(
 153       [](auto& a, const auto& b) { a = b; },
 154       [](auto&, const auto&) {
 155         BOOST_THROW_EXCEPTION(
 156             std::invalid_argument("cannot assign axes, types do not match"));
 157       },
 158       t, u);
 159 }
 160
 161 template <class... Ts, class U>
 162 void axes_assign(std::tuple<Ts...>& t, const U& u) {
 163   using namespace ::boost::mp11;
 164   mp_for_each<mp_iota_c<sizeof...(Ts)>>([&](auto I) {
 165     using T = mp_at_c<std::tuple<Ts...>, I>;
 166     std::get<I>(t) = axis::get<T>(u[I]);
 167   });
 168 }
 169
 170 template <class T, class... Us>
 171 void axes_assign(T& t, const std::tuple<Us...>& u) {
 172   // resize instead of reserve, because t may not be empty and we want exact capacity
 173   t.resize(sizeof...(Us));
 174   using namespace ::boost::mp11;
 175   mp_for_each<mp_iota_c<sizeof...(Us)>>([&](auto I) { t[I] = std::get<I>(u); });
 176 }
 177
 178 template <typename T, typename U>
 179 void axes_assign(T& t, const U& u) {
 180   t.assign(u.begin(), u.end());
 181 }
 182
 183 template <class Archive, class T>
 184 void axes_serialize(Archive& ar, T& axes) {
 185   ar& make_nvp("axes", axes);
 186 }
 187
 188 template <class Archive, class... Ts>
 189 void axes_serialize(Archive& ar, std::tuple<Ts...>& axes) {
 190   // needed to keep serialization format backward compatible
 191   struct proxy {
 192     std::tuple<Ts...>& t;
 193     void serialize(Archive& ar, unsigned /* version */) {
 194       mp11::tuple_for_each(t, [&ar](auto& x) { ar& make_nvp("item", x); });
 195     }
 196   };
 197   proxy p{axes};
 198   ar& make_nvp("axes", p);
 199 }
 200
 201 // create empty dynamic axis which can store any axes types from the argument
 202 template <class T>
 203 auto make_empty_dynamic_axes(const T& axes) {
 204   return make_default(axes);
 205 }
 206
 207 template <class... Ts>
 208 auto make_empty_dynamic_axes(const std::tuple<Ts...>&) {
 209   using namespace ::boost::mp11;
 210   using L = mp_unique<axis::variant<Ts...>>;
 211   // return std::vector<axis::variant<Axis0, Axis1, ...>> or std::vector<Axis0>
 212   return std::vector<mp_if_c<(mp_size<L>::value == 1), mp_first<L>, L>>{};
 213 }
 214
 215 template <class T>
 216 void axis_index_is_valid(const T& axes, const unsigned N) {
 217   BOOST_ASSERT_MSG(N < axes_rank(axes), "index out of range");
 218 }
 219
 220 template <class Axes, class V>
 221 void for_each_axis_impl(std::true_type, Axes&& axes, V&& v) {
 222   for (auto&& a : axes) { axis::visit(std::forward<V>(v), a); }
 223 }
 224
 225 template <class Axes, class V>
 226 void for_each_axis_impl(std::false_type, Axes&& axes, V&& v) {
 227   for (auto&& a : axes) std::forward<V>(v)(a);
 228 }
 229
 230 template <class Axes, class V>
 231 void for_each_axis(Axes&& a, V&& v) {
 232   using namespace ::boost::mp11;
 233   using T = mp_first<std::decay_t<Axes>>;
 234   for_each_axis_impl(is_axis_variant<T>(), std::forward<Axes>(a), std::forward<V>(v));
 235 }
 236
 237 template <class V, class... Axis>
 238 void for_each_axis(const std::tuple<Axis...>& a, V&& v) {
 239   mp11::tuple_for_each(a, std::forward<V>(v));
 240 }
 241
 242 template <class V, class... Axis>
 243 void for_each_axis(std::tuple<Axis...>& a, V&& v) {
 244   mp11::tuple_for_each(a, std::forward<V>(v));
 245 }
 246
 247 // total number of bins including *flow bins
 248 template <class T>
 249 std::size_t bincount(const T& axes) {
 250   std::size_t n = 1;
 251   for_each_axis(axes, [&n](const auto& a) {
 252     const auto old = n;
 253     const auto s = axis::traits::extent(a);
 254     n *= s;
 255     if (s > 0 && n < old) BOOST_THROW_EXCEPTION(std::overflow_error("bincount overflow"));
 256   });
 257   return n;
 258 }
 259
 260 // initial offset for the linear index
 261 template <class T>
 262 std::size_t offset(const T& axes) {
 263   std::size_t n = 0;
 264   for_each_axis(axes, [&n, stride = static_cast<std::size_t>(1)](const auto& a) mutable {
 265     if (axis::traits::options(a) & axis::option::growth)
 266       n = invalid_index;
 267     else if (n != invalid_index && axis::traits::options(a) & axis::option::underflow)
 268       n += stride;
 269     stride *= axis::traits::extent(a);
 270   });
 271   return n;
 272 }
 273
 274 template <class T>
 275 using buffer_size_impl = typename std::tuple_size<T>::type;
 276
 277 template <class T>
 278 using buffer_size = mp11::mp_eval_or<
 279     std::integral_constant<std::size_t, BOOST_HISTOGRAM_DETAIL_AXES_LIMIT>,
 280     buffer_size_impl, T>;
 281
 282 template <class T, std::size_t N>
 283 class sub_array : public std::array<T, N> {
 284   using base_type = std::array<T, N>;
 285
 286 public:
 287   explicit sub_array(std::size_t s) noexcept(
 288       std::is_nothrow_default_constructible<T>::value)
 289       : size_(s) {
 290     BOOST_ASSERT_MSG(size_ <= N, "requested size exceeds size of static buffer");
 291   }
 292
 293   sub_array(std::size_t s,
 294             const T& value) noexcept(std::is_nothrow_copy_constructible<T>::value)
 295       : size_(s) {
 296     BOOST_ASSERT_MSG(size_ <= N, "requested size exceeds size of static buffer");
 297     std::array<T, N>::fill(value);
 298   }
 299
 300   // need to override both versions of std::array
 301   auto end() noexcept { return base_type::begin() + size_; }
 302   auto end() const noexcept { return base_type::begin() + size_; }
 303
 304   auto size() const noexcept { return size_; }
 305
 306 private:
 307   std::size_t size_;
 308 };
 309
 310 template <class U, class T>
 311 using stack_buffer = sub_array<U, buffer_size<T>::value>;
 312
 313 // make default-constructed buffer (no initialization for POD types)
 314 template <class U, class T>
 315 auto make_stack_buffer(const T& t) {
 316   return stack_buffer<U, T>(axes_rank(t));
 317 }
 318
 319 // make buffer with elements initialized to v
 320 template <class U, class T, class V>
 321 auto make_stack_buffer(const T& t, V&& v) {
 322   return stack_buffer<U, T>(axes_rank(t), std::forward<V>(v));
 323 }
 324
 325 template <class T>
 326 using has_underflow =
 327     decltype(axis::traits::static_options<T>::test(axis::option::underflow));
 328
 329 template <class T>
 330 using is_growing = decltype(axis::traits::static_options<T>::test(axis::option::growth));
 331
 332 template <class T>
 333 using is_not_inclusive = mp11::mp_not<axis::traits::static_is_inclusive<T>>;
 334
 335 // for vector<T>
 336 template <class T>
 337 struct axis_types_impl {
 338   using type = mp11::mp_list<std::decay_t<T>>;
 339 };
 340
 341 // for vector<variant<Ts...>>
 342 template <class... Ts>
 343 struct axis_types_impl<axis::variant<Ts...>> {
 344   using type = mp11::mp_list<std::decay_t<Ts>...>;
 345 };
 346
 347 // for tuple<Ts...>
 348 template <class... Ts>
 349 struct axis_types_impl<std::tuple<Ts...>> {
 350   using type = mp11::mp_list<std::decay_t<Ts>...>;
 351 };
 352
 353 template <class T>
 354 using axis_types =
 355     typename axis_types_impl<mp11::mp_if<is_vector_like<T>, mp11::mp_first<T>, T>>::type;
 356
 357 template <template <class> class Trait, class Axes>
 358 using has_special_axis = mp11::mp_any_of<axis_types<Axes>, Trait>;
 359
 360 template <class Axes>
 361 using has_growing_axis = mp11::mp_any_of<axis_types<Axes>, is_growing>;
 362
 363 template <class Axes>
 364 using has_non_inclusive_axis = mp11::mp_any_of<axis_types<Axes>, is_not_inclusive>;
 365
 366 template <class T>
 367 constexpr std::size_t type_score() {
 368   return sizeof(T) *
 369          (std::is_integral<T>::value ? 1 : std::is_floating_point<T>::value ? 10 : 100);
 370 }
 371
 372 // arbitrary ordering of types
 373 template <class T, class U>
 374 using type_less = mp11::mp_bool<(type_score<T>() < type_score<U>())>;
 375
 376 template <class Axes>
 377 using value_types = mp11::mp_sort<
 378     mp11::mp_unique<mp11::mp_transform<axis::traits::value_type, axis_types<Axes>>>,
 379     type_less>;
 380
 381 } // namespace detail
 382 } // namespace histogram
 383 } // namespace boost
 384
 385 #endif