project.hpp 3.6 KB

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