histogram/detail/axes.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_AXES_HPP
0008 #define BOOST_HISTOGRAM_DETAIL_AXES_HPP
0009
0010 #include <array>
0011 #include <boost/core/nvp.hpp>
0012 #include <boost/histogram/axis/traits.hpp>
0013 #include <boost/histogram/axis/variant.hpp>
0014 #include <boost/histogram/detail/make_default.hpp>
0015 #include <boost/histogram/detail/nonmember_container_access.hpp>
0016 #include <boost/histogram/detail/optional_index.hpp>
0017 #include <boost/histogram/detail/priority.hpp>
0018 #include <boost/histogram/detail/relaxed_tuple_size.hpp>
0019 #include <boost/histogram/detail/static_if.hpp>
0020 #include <boost/histogram/detail/sub_array.hpp>
0021 #include <boost/histogram/detail/try_cast.hpp>
0022 #include <boost/histogram/fwd.hpp>
0023 #include <boost/mp11/algorithm.hpp>
0024 #include <boost/mp11/integer_sequence.hpp>
0025 #include <boost/mp11/list.hpp>
0026 #include <boost/mp11/tuple.hpp>
0027 #include <boost/mp11/utility.hpp>
0028 #include <boost/throw_exception.hpp>
0029 #include <cassert>
0030 #include <initializer_list>
0031 #include <iterator>
0032 #include <stdexcept>
0033 #include <string>
0034 #include <tuple>
0035 #include <type_traits>
0036 #include <vector>
0037
0038 namespace boost {
0039 namespace histogram {
0040 namespace detail {
0041
0042 template <class T, class Unary>
0043 void for_each_axis_impl(dynamic_size, T& t, Unary& p) {
0044   for (auto& a : t) axis::visit(p, a);
0045 }
0046
0047 template <class N, class T, class Unary>
0048 void for_each_axis_impl(N, T& t, Unary& p) {
0049   mp11::tuple_for_each(t, p);
0050 }
0051
0052 // also matches const T and const Unary
0053 template <class T, class Unary>
0054 void for_each_axis(T&& t, Unary&& p) {
0055   for_each_axis_impl(relaxed_tuple_size(t), t, p);
0056 }
0057
0058 // merge if a and b are discrete and growing
0059 struct axis_merger {
0060   template <class T, class U>
0061   T operator()(const T& a, const U& u) {
0062     const T* bp = ptr_cast<T>(&u);
0063     if (!bp) BOOST_THROW_EXCEPTION(std::invalid_argument("axes not mergable"));
0064     using O = axis::traits::get_options<T>;
0065     constexpr bool discrete_and_growing =
0066         axis::traits::is_continuous<T>::value == false && O::test(axis::option::growth);
0067     return impl(mp11::mp_bool<discrete_and_growing>{}, a, *bp);
0068   }
0069
0070   template <class T>
0071   T impl(std::false_type, const T& a, const T& b) {
0072     if (!relaxed_equal{}(a, b))
0073       BOOST_THROW_EXCEPTION(std::invalid_argument("axes not mergable"));
0074     return a;
0075   }
0076
0077   template <class T>
0078   T impl(std::true_type, const T& a, const T& b) {
0079     if (relaxed_equal{}(axis::traits::metadata(a), axis::traits::metadata(b))) {
0080       auto r = a;
0081       if (axis::traits::is_ordered<T>::value) {
0082         r.update(b.value(0));
0083         r.update(b.value(b.size() - 1));
0084       } else
0085         for (auto&& v : b) r.update(v);
0086       return r;
0087     }
0088     return impl(std::false_type{}, a, b);
0089   }
0090 };
0091
0092 // create empty dynamic axis which can store any axes types from the argument
0093 template <class T>
0094 auto make_empty_dynamic_axes(const T& axes) {
0095   return make_default(axes);
0096 }
0097
0098 template <class... Ts>
0099 auto make_empty_dynamic_axes(const std::tuple<Ts...>&) {
0100   using namespace ::boost::mp11;
0101   using L = mp_unique<axis::variant<Ts...>>;
0102   // return std::vector<axis::variant<Axis0, Axis1, ...>> or std::vector<Axis0>
0103   return std::vector<mp_if_c<(mp_size<L>::value == 1), mp_first<L>, L>>{};
0104 }
0105
0106 template <class T, class Functor, std::size_t... Is>
0107 auto axes_transform_impl(const T& t, Functor&& f, mp11::index_sequence<Is...>) {
0108   return std::make_tuple(f(Is, std::get<Is>(t))...);
0109 }
0110
0111 // warning: sequential order of functor execution is platform-dependent!
0112 template <class... Ts, class Functor>
0113 auto axes_transform(const std::tuple<Ts...>& old_axes, Functor&& f) {
0114   return axes_transform_impl(old_axes, std::forward<Functor>(f),
0115                              mp11::make_index_sequence<sizeof...(Ts)>{});
0116 }
0117
0118 // changing axes type is not supported
0119 template <class T, class Functor>
0120 T axes_transform(const T& old_axes, Functor&& f) {
0121   T axes = make_default(old_axes);
0122   axes.reserve(old_axes.size());
0123   for_each_axis(old_axes, [&](const auto& a) { axes.emplace_back(f(axes.size(), a)); });
0124   return axes;
0125 }
0126
0127 template <class... Ts, class Binary, std::size_t... Is>
0128 std::tuple<Ts...> axes_transform_impl(const std::tuple<Ts...>& lhs,
0129                                       const std::tuple<Ts...>& rhs, Binary&& bin,
0130                                       mp11::index_sequence<Is...>) {
0131   return std::make_tuple(bin(std::get<Is>(lhs), std::get<Is>(rhs))...);
0132 }
0133
0134 template <class... Ts, class Binary>
0135 std::tuple<Ts...> axes_transform(const std::tuple<Ts...>& lhs,
0136                                  const std::tuple<Ts...>& rhs, Binary&& bin) {
0137   return axes_transform_impl(lhs, rhs, bin, mp11::make_index_sequence<sizeof...(Ts)>{});
0138 }
0139
0140 template <class T, class Binary>
0141 T axes_transform(const T& lhs, const T& rhs, Binary&& bin) {
0142   T ax = make_default(lhs);
0143   ax.reserve(lhs.size());
0144   using std::begin;
0145   auto ir = begin(rhs);
0146   for (auto&& li : lhs) {
0147     axis::visit(
0148         [&](const auto& li) {
0149           axis::visit([&](const auto& ri) { ax.emplace_back(bin(li, ri)); }, *ir);
0150         },
0151         li);
0152     ++ir;
0153   }
0154   return ax;
0155 }
0156
0157 template <class T>
0158 unsigned axes_rank(const T& axes) {
0159   using std::begin;
0160   using std::end;
0161   return static_cast<unsigned>(std::distance(begin(axes), end(axes)));
0162 }
0163
0164 template <class... Ts>
0165 constexpr unsigned axes_rank(const std::tuple<Ts...>&) {
0166   return static_cast<unsigned>(sizeof...(Ts));
0167 }
0168
0169 template <class T>
0170 void throw_if_axes_is_too_large(const T& axes) {
0171   if (axes_rank(axes) > BOOST_HISTOGRAM_DETAIL_AXES_LIMIT)
0172     BOOST_THROW_EXCEPTION(
0173         std::invalid_argument("length of axis vector exceeds internal buffers, "
0174                               "recompile with "
0175                               "-DBOOST_HISTOGRAM_DETAIL_AXES_LIMIT=<new max size> "
0176                               "to increase internal buffers"));
0177 }
0178
0179 // tuple is never too large because internal buffers adapt to size of tuple
0180 template <class... Ts>
0181 void throw_if_axes_is_too_large(const std::tuple<Ts...>&) {}
0182
0183 template <unsigned N, class... Ts>
0184 decltype(auto) axis_get(std::tuple<Ts...>& axes) {
0185   return std::get<N>(axes);
0186 }
0187
0188 template <unsigned N, class... Ts>
0189 decltype(auto) axis_get(const std::tuple<Ts...>& axes) {
0190   return std::get<N>(axes);
0191 }
0192
0193 template <unsigned N, class T>
0194 decltype(auto) axis_get(T& axes) {
0195   return axes[N];
0196 }
0197
0198 template <unsigned N, class T>
0199 decltype(auto) axis_get(const T& axes) {
0200   return axes[N];
0201 }
0202
0203 template <class... Ts>
0204 auto axis_get(std::tuple<Ts...>& axes, const unsigned i) {
0205   constexpr auto S = sizeof...(Ts);
0206   using V = mp11::mp_unique<axis::variant<Ts*...>>;
0207   return mp11::mp_with_index<S>(i, [&axes](auto i) { return V(&std::get<i>(axes)); });
0208 }
0209
0210 template <class... Ts>
0211 auto axis_get(const std::tuple<Ts...>& axes, const unsigned i) {
0212   constexpr auto S = sizeof...(Ts);
0213   using V = mp11::mp_unique<axis::variant<const Ts*...>>;
0214   return mp11::mp_with_index<S>(i, [&axes](auto i) { return V(&std::get<i>(axes)); });
0215 }
0216
0217 template <class T>
0218 decltype(auto) axis_get(T& axes, const unsigned i) {
0219   return axes[i];
0220 }
0221
0222 template <class T>
0223 decltype(auto) axis_get(const T& axes, const unsigned i) {
0224   return axes[i];
0225 }
0226
0227 template <class T, class U, std::size_t... Is>
0228 bool axes_equal_impl(const T& t, const U& u, mp11::index_sequence<Is...>) noexcept {
0229   bool result = true;
0230   // operator folding emulation
0231   (void)std::initializer_list<bool>{
0232       (result &= relaxed_equal{}(std::get<Is>(t), std::get<Is>(u)))...};
0233   return result;
0234 }
0235
0236 template <class... Ts, class... Us>
0237 bool axes_equal_impl(const std::tuple<Ts...>& t, const std::tuple<Us...>& u) noexcept {
0238   return axes_equal_impl(
0239       t, u, mp11::make_index_sequence<(std::min)(sizeof...(Ts), sizeof...(Us))>{});
0240 }
0241
0242 template <class... Ts, class U>
0243 bool axes_equal_impl(const std::tuple<Ts...>& t, const U& u) noexcept {
0244   using std::begin;
0245   auto iu = begin(u);
0246   bool result = true;
0247   mp11::tuple_for_each(t, [&](const auto& ti) {
0248     axis::visit([&](const auto& ui) { result &= relaxed_equal{}(ti, ui); }, *iu);
0249     ++iu;
0250   });
0251   return result;
0252 }
0253
0254 template <class T, class... Us>
0255 bool axes_equal_impl(const T& t, const std::tuple<Us...>& u) noexcept {
0256   return axes_equal_impl(u, t);
0257 }
0258
0259 template <class T, class U>
0260 bool axes_equal_impl(const T& t, const U& u) noexcept {
0261   using std::begin;
0262   auto iu = begin(u);
0263   bool result = true;
0264   for (auto&& ti : t) {
0265     axis::visit(
0266         [&](const auto& ti) {
0267           axis::visit([&](const auto& ui) { result &= relaxed_equal{}(ti, ui); }, *iu);
0268         },
0269         ti);
0270     ++iu;
0271   }
0272   return result;
0273 }
0274
0275 template <class T, class U>
0276 bool axes_equal(const T& t, const U& u) noexcept {
0277   return axes_rank(t) == axes_rank(u) && axes_equal_impl(t, u);
0278 }
0279
0280 // enable_if_t needed by msvc :(
0281 template <class... Ts, class... Us>
0282 std::enable_if_t<!(std::is_same<std::tuple<Ts...>, std::tuple<Us...>>::value)>
0283 axes_assign(std::tuple<Ts...>&, const std::tuple<Us...>&) {
0284   BOOST_THROW_EXCEPTION(std::invalid_argument("cannot assign axes, types do not match"));
0285 }
0286
0287 template <class... Ts>
0288 void axes_assign(std::tuple<Ts...>& t, const std::tuple<Ts...>& u) {
0289   t = u;
0290 }
0291
0292 template <class... Ts, class U>
0293 void axes_assign(std::tuple<Ts...>& t, const U& u) {
0294   if (sizeof...(Ts) == detail::size(u)) {
0295     using std::begin;
0296     auto iu = begin(u);
0297     mp11::tuple_for_each(t, [&](auto& ti) {
0298       using T = std::decay_t<decltype(ti)>;
0299       ti = axis::get<T>(*iu);
0300       ++iu;
0301     });
0302     return;
0303   }
0304   BOOST_THROW_EXCEPTION(std::invalid_argument("cannot assign axes, sizes do not match"));
0305 }
0306
0307 template <class T, class... Us>
0308 void axes_assign(T& t, const std::tuple<Us...>& u) {
0309   // resize instead of reserve, because t may not be empty and we want exact capacity
0310   t.resize(sizeof...(Us));
0311   using std::begin;
0312   auto it = begin(t);
0313   mp11::tuple_for_each(u, [&](const auto& ui) { *it++ = ui; });
0314 }
0315
0316 template <class T, class U>
0317 void axes_assign(T& t, const U& u) {
0318   t.assign(u.begin(), u.end());
0319 }
0320
0321 template <class Archive, class T>
0322 void axes_serialize(Archive& ar, T& axes) {
0323   ar& make_nvp("axes", axes);
0324 }
0325
0326 template <class Archive, class... Ts>
0327 void axes_serialize(Archive& ar, std::tuple<Ts...>& axes) {
0328   // needed to keep serialization format backward compatible
0329   struct proxy {
0330     std::tuple<Ts...>& t;
0331     void serialize(Archive& ar, unsigned /* version */) {
0332       mp11::tuple_for_each(t, [&ar](auto& x) { ar& make_nvp("item", x); });
0333     }
0334   };
0335   proxy p{axes};
0336   ar& make_nvp("axes", p);
0337 }
0338
0339 // total number of bins including *flow bins
0340 template <class T>
0341 std::size_t bincount(const T& axes) {
0342   std::size_t n = 1;
0343   for_each_axis(axes, [&n](const auto& a) {
0344     const auto old = n;
0345     const auto s = axis::traits::extent(a);
0346     n *= s;
0347     if (s > 0 && n < old) BOOST_THROW_EXCEPTION(std::overflow_error("bincount overflow"));
0348   });
0349   return n;
0350 }
0351
0352 /** Initial offset for the linear index.
0353
0354   This precomputes an offset for the global index so that axis index = -1 addresses the
0355   first entry in the storage. Example: one-dim. histogram. The offset is 1, stride is 1,
0356   and global_index = offset + axis_index * stride == 0 addresses the first element of
0357   the storage.
0358
0359   Using the offset makes the hot inner loop that computes the global_index simpler and
0360   thus faster, because we do not have to branch for each axis to check whether it has
0361   an underflow bin.
0362
0363   The offset is set to an invalid value when the histogram contains at least one growing
0364   axis, because this optimization then cannot be used. See detail/linearize.hpp, in this
0365   case linearize_growth is called.
0366 */
0367 template <class T>
0368 std::size_t offset(const T& axes) {
0369   std::size_t n = 0;
0370   auto stride = static_cast<std::size_t>(1);
0371   for_each_axis(axes, [&](const auto& a) {
0372     if (axis::traits::options(a) & axis::option::growth)
0373       n = invalid_index;
0374     else if (n != invalid_index && axis::traits::options(a) & axis::option::underflow)
0375       n += stride;
0376     stride *= axis::traits::extent(a);
0377   });
0378   return n;
0379 }
0380
0381 // make default-constructed buffer (no initialization for POD types)
0382 template <class T, class A>
0383 auto make_stack_buffer(const A& a) {
0384   return sub_array<T, buffer_size<A>::value>(axes_rank(a));
0385 }
0386
0387 // make buffer with elements initialized to v
0388 template <class T, class A>
0389 auto make_stack_buffer(const A& a, const T& t) {
0390   return sub_array<T, buffer_size<A>::value>(axes_rank(a), t);
0391 }
0392
0393 template <class T>
0394 using has_underflow =
0395     decltype(axis::traits::get_options<T>::test(axis::option::underflow));
0396
0397 template <class T>
0398 using is_growing = decltype(axis::traits::get_options<T>::test(axis::option::growth));
0399
0400 template <class T>
0401 using is_not_inclusive = mp11::mp_not<axis::traits::is_inclusive<T>>;
0402
0403 // for vector<T>
0404 template <class T>
0405 struct axis_types_impl {
0406   using type = mp11::mp_list<std::decay_t<T>>;
0407 };
0408
0409 // for vector<variant<Ts...>>
0410 template <class... Ts>
0411 struct axis_types_impl<axis::variant<Ts...>> {
0412   using type = mp11::mp_list<std::decay_t<Ts>...>;
0413 };
0414
0415 // for tuple<Ts...>
0416 template <class... Ts>
0417 struct axis_types_impl<std::tuple<Ts...>> {
0418   using type = mp11::mp_list<std::decay_t<Ts>...>;
0419 };
0420
0421 template <class T>
0422 using axis_types =
0423     typename axis_types_impl<mp11::mp_if<is_vector_like<T>, mp11::mp_first<T>, T>>::type;
0424
0425 template <template <class> class Trait, class Axes>
0426 using has_special_axis = mp11::mp_any_of<axis_types<Axes>, Trait>;
0427
0428 template <class Axes>
0429 using has_growing_axis = mp11::mp_any_of<axis_types<Axes>, is_growing>;
0430
0431 template <class Axes>
0432 using has_non_inclusive_axis = mp11::mp_any_of<axis_types<Axes>, is_not_inclusive>;
0433
0434 template <class T>
0435 constexpr std::size_t type_score() {
0436   return sizeof(T) * (std::is_integral<T>::value         ? 1
0437                       : std::is_floating_point<T>::value ? 10
0438                                                          : 100);
0439 }
0440
0441 // arbitrary ordering of types
0442 template <class T, class U>
0443 using type_less = mp11::mp_bool<(type_score<T>() < type_score<U>())>;
0444
0445 template <class Axes>
0446 using value_types = mp11::mp_sort<
0447     mp11::mp_unique<mp11::mp_transform<axis::traits::value_type, axis_types<Axes>>>,
0448     type_less>;
0449
0450 } // namespace detail
0451 } // namespace histogram
0452 } // namespace boost
0453
0454 #endif