boost/histogram/algorithm/project.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_ALGORITHM_PROJECT_HPP
   8 #define BOOST_HISTOGRAM_ALGORITHM_PROJECT_HPP
   9
  10 #include <algorithm>
  11 #include <boost/histogram/axis/variant.hpp>
  12 #include <boost/histogram/detail/detect.hpp>
  13 #include <boost/histogram/detail/make_default.hpp>
  14 #include <boost/histogram/detail/static_if.hpp>
  15 #include <boost/histogram/histogram.hpp>
  16 #include <boost/histogram/indexed.hpp>
  17 #include <boost/histogram/unsafe_access.hpp>
  18 #include <boost/mp11/list.hpp>
  19 #include <boost/mp11/set.hpp>
  20 #include <boost/mp11/utility.hpp>
  21 #include <boost/throw_exception.hpp>
  22 #include <stdexcept>
  23 #include <type_traits>
  24 #include <vector>
  25
  26 namespace boost {
  27 namespace histogram {
  28 namespace algorithm {
  29
  30 /**
  31   Returns a lower-dimensional histogram, summing over removed axes.
  32
  33   Arguments are the source histogram and compile-time numbers, the remaining indices of
  34   the axes. Returns a new histogram which only contains the subset of axes. The source
  35   histogram is summed over the removed axes.
  36 */
  37 template <class A, class S, unsigned N, typename... Ns>
  38 auto project(const histogram<A, S>& h, std::integral_constant<unsigned, N>, Ns...) {
  39   using LN = mp11::mp_list<std::integral_constant<unsigned, N>, Ns...>;
  40   static_assert(mp11::mp_is_set<LN>::value, "indices must be unique");
  41
  42   const auto& old_axes = unsafe_access::axes(h);
  43   auto axes = detail::static_if<detail::is_tuple<A>>(
  44       [&](const auto& old_axes) {
  45         return std::make_tuple(std::get<N>(old_axes), std::get<Ns::value>(old_axes)...);
  46       },
  47       [&](const auto& old_axes) {
  48         return std::decay_t<decltype(old_axes)>({old_axes[N], old_axes[Ns::value]...});
  49       },
  50       old_axes);
  51
  52   const auto& old_storage = unsafe_access::storage(h);
  53   using A2 = decltype(axes);
  54   auto result = histogram<A2, S>(std::move(axes), detail::make_default(old_storage));
  55   auto idx = detail::make_stack_buffer<int>(unsafe_access::axes(result));
  56   for (auto&& x : indexed(h, coverage::all)) {
  57     auto i = idx.begin();
  58     mp11::mp_for_each<LN>([&i, &x](auto J) { *i++ = x.index(J); });
  59     result.at(idx) += *x;
  60   }
  61   return result;
  62 }
  63
  64 /**
  65   Returns a lower-dimensional histogram, summing over removed axes.
  66
  67   This version accepts a source histogram and an iterable range containing the remaining
  68   indices.
  69 */
  70 template <class A, class S, class Iterable, class = detail::requires_iterable<Iterable>>
  71 auto project(const histogram<A, S>& h, const Iterable& c) {
  72   using namespace boost::mp11;
  73   const auto& old_axes = unsafe_access::axes(h);
  74
  75   // axes is always std::vector<...>, even if A is tuple
  76   auto axes = detail::make_empty_dynamic_axes(old_axes);
  77   axes.reserve(c.size());
  78   auto seen = detail::make_stack_buffer<bool>(old_axes, false);
  79   for (auto d : c) {
  80     if (static_cast<unsigned>(d) >= h.rank())
  81       BOOST_THROW_EXCEPTION(std::invalid_argument("invalid axis index"));
  82     if (seen[d]) BOOST_THROW_EXCEPTION(std::invalid_argument("indices are not unique"));
  83     seen[d] = true;
  84     axes.emplace_back(detail::axis_get(old_axes, d));
  85   }
  86
  87   const auto& old_storage = unsafe_access::storage(h);
  88   auto result =
  89       histogram<decltype(axes), S>(std::move(axes), detail::make_default(old_storage));
  90   auto idx = detail::make_stack_buffer<int>(unsafe_access::axes(result));
  91   for (auto&& x : indexed(h, coverage::all)) {
  92     auto i = idx.begin();
  93     for (auto d : c) *i++ = x.index(d);
  94     result.at(idx) += *x;
  95   }
  96
  97   return result;
  98 }
  99
 100 } // namespace algorithm
 101 } // namespace histogram
 102 } // namespace boost
 103
 104 #endif