histogram/algorithm/project.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_ALGORITHM_PROJECT_HPP
0008 #define BOOST_HISTOGRAM_ALGORITHM_PROJECT_HPP
0009
0010 #include <algorithm>
0011 #include <boost/histogram/axis/variant.hpp>
0012 #include <boost/histogram/detail/detect.hpp>
0013 #include <boost/histogram/detail/make_default.hpp>
0014 #include <boost/histogram/detail/static_if.hpp>
0015 #include <boost/histogram/histogram.hpp>
0016 #include <boost/histogram/indexed.hpp>
0017 #include <boost/histogram/unsafe_access.hpp>
0018 #include <boost/mp11/list.hpp>
0019 #include <boost/mp11/set.hpp>
0020 #include <boost/mp11/utility.hpp>
0021 #include <boost/throw_exception.hpp>
0022 #include <stdexcept>
0023 #include <type_traits>
0024 #include <vector>
0025
0026 namespace boost {
0027 namespace histogram {
0028 namespace algorithm {
0029
0030 /**
0031   Returns a lower-dimensional histogram, summing over removed axes.
0032
0033   Arguments are the source histogram and compile-time numbers, the remaining indices of
0034   the axes. Returns a new histogram which only contains the subset of axes. The source
0035   histogram is summed over the removed axes.
0036 */
0037 template <class A, class S, unsigned N, typename... Ns>
0038 auto project(const histogram<A, S>& h, std::integral_constant<unsigned, N>, Ns...) {
0039   using LN = mp11::mp_list<std::integral_constant<unsigned, N>, Ns...>;
0040   static_assert(mp11::mp_is_set<LN>::value, "indices must be unique");
0041
0042   const auto& old_axes = unsafe_access::axes(h);
0043   auto axes = detail::static_if<detail::is_tuple<A>>(
0044       [&](const auto& old_axes) {
0045         return std::make_tuple(std::get<N>(old_axes), std::get<Ns::value>(old_axes)...);
0046       },
0047       [&](const auto& old_axes) {
0048         return std::decay_t<decltype(old_axes)>({old_axes[N], old_axes[Ns::value]...});
0049       },
0050       old_axes);
0051
0052   const auto& old_storage = unsafe_access::storage(h);
0053   using A2 = decltype(axes);
0054   auto result = histogram<A2, S>(std::move(axes), detail::make_default(old_storage));
0055   auto idx = detail::make_stack_buffer<int>(unsafe_access::axes(result));
0056   for (auto&& x : indexed(h, coverage::all)) {
0057     auto i = idx.begin();
0058     mp11::mp_for_each<LN>([&i, &x](auto J) { *i++ = x.index(J); });
0059     result.at(idx) += *x;
0060   }
0061   return result;
0062 }
0063
0064 /**
0065   Returns a lower-dimensional histogram, summing over removed axes.
0066
0067   This version accepts a source histogram and an iterable range containing the remaining
0068   indices.
0069 */
0070 template <class A, class S, class Iterable, class = detail::requires_iterable<Iterable>>
0071 auto project(const histogram<A, S>& h, const Iterable& c) {
0072   using namespace boost::mp11;
0073   const auto& old_axes = unsafe_access::axes(h);
0074
0075   // axes is always std::vector<...>, even if A is tuple
0076   auto axes = detail::make_empty_dynamic_axes(old_axes);
0077   axes.reserve(c.size());
0078   auto seen = detail::make_stack_buffer(old_axes, false);
0079   for (auto d : c) {
0080     if (static_cast<unsigned>(d) >= h.rank())
0081       BOOST_THROW_EXCEPTION(std::invalid_argument("invalid axis index"));
0082     if (seen[d]) BOOST_THROW_EXCEPTION(std::invalid_argument("indices are not unique"));
0083     seen[d] = true;
0084     axes.emplace_back(detail::axis_get(old_axes, d));
0085   }
0086
0087   const auto& old_storage = unsafe_access::storage(h);
0088   auto result =
0089       histogram<decltype(axes), S>(std::move(axes), detail::make_default(old_storage));
0090   auto idx = detail::make_stack_buffer<int>(unsafe_access::axes(result));
0091   for (auto&& x : indexed(h, coverage::all)) {
0092     auto i = idx.begin();
0093     for (auto d : c) *i++ = x.index(d);
0094     result.at(idx) += *x;
0095   }
0096
0097   return result;
0098 }
0099
0100 } // namespace algorithm
0101 } // namespace histogram
0102 } // namespace boost
0103
0104 #endif