reduce.hpp 6.7 KB

123456789101112131415161718192021222324252627282930313233343536373839404142434445464748495051525354555657585960616263646566676869707172737475767778798081828384858687888990919293949596979899100101102103104105106107108109110111112113114115116117118119120121122123124125126127128129130131132133134135136137138139140141142143144145146147148149150151152153154155156157158159160161162163164165166167168169170171172173174175176177178179180181182183184185186187188189190191192193194195196197198199200201202203204205206207208209210211212213214215216217218219220221222223
  1. // Copyright 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_REDUCE_HPP
  7. #define BOOST_HISTOGRAM_ALGORITHM_REDUCE_HPP
  8. #include <boost/assert.hpp>
  9. #include <boost/histogram/detail/axes.hpp>
  10. #include <boost/histogram/detail/meta.hpp>
  11. #include <boost/histogram/fwd.hpp>
  12. #include <boost/histogram/indexed.hpp>
  13. #include <boost/histogram/unsafe_access.hpp>
  14. #include <boost/throw_exception.hpp>
  15. #include <cmath>
  16. #include <limits>
  17. #include <stdexcept>
  18. #include <type_traits>
  19. namespace boost {
  20. namespace histogram {
  21. namespace algorithm {
  22. /**
  23. Option type returned by the helper functions shrink_and_rebin(), shrink(), rebin().
  24. */
  25. struct reduce_option {
  26. #ifndef BOOST_HISTOGRAM_DOXYGEN_INVOKED
  27. unsigned iaxis = 0;
  28. double lower, upper;
  29. unsigned merge = 0;
  30. reduce_option() noexcept = default;
  31. reduce_option(unsigned i, double l, double u, unsigned m)
  32. : iaxis(i), lower(l), upper(u), merge(m) {
  33. if (lower == upper)
  34. BOOST_THROW_EXCEPTION(std::invalid_argument("lower != upper required"));
  35. if (merge == 0) BOOST_THROW_EXCEPTION(std::invalid_argument("merge > 0 required"));
  36. }
  37. #endif
  38. };
  39. /**
  40. Shrink and rebin option.
  41. @param iaxis which axis to operate on.
  42. @param lower lowest bound that should be kept.
  43. @param upper highest bound that should be kept. If upper is inside bin interval, the
  44. whole interval is removed.
  45. @param merge how many adjacent bins to merge into one.
  46. */
  47. inline auto shrink_and_rebin(unsigned iaxis, double lower, double upper, unsigned merge) {
  48. return reduce_option{iaxis, lower, upper, merge};
  49. }
  50. /**
  51. Shrink option.
  52. @param iaxis which axis to operate on.
  53. @param lower lowest bound that should be kept.
  54. @param upper highest bound that should be kept. If upper is inside bin interval, the
  55. whole interval is removed.
  56. */
  57. inline auto shrink(unsigned iaxis, double lower, double upper) {
  58. return reduce_option{iaxis, lower, upper, 1};
  59. }
  60. /**
  61. Rebin option.
  62. @param iaxis which axis to operate on.
  63. @param merge how many adjacent bins to merge into one.
  64. */
  65. inline auto rebin(unsigned iaxis, unsigned merge) {
  66. return reduce_option{iaxis, std::numeric_limits<double>::quiet_NaN(),
  67. std::numeric_limits<double>::quiet_NaN(), merge};
  68. }
  69. /**
  70. Convenience overload for single axis.
  71. @param lower lowest bound that should be kept.
  72. @param upper highest bound that should be kept. If upper is inside bin interval, the
  73. whole interval is removed.
  74. @param merge how many adjacent bins to merge into one.
  75. */
  76. inline auto shrink_and_rebin(double lower, double upper, unsigned merge) {
  77. return shrink_and_rebin(0, lower, upper, merge);
  78. }
  79. /**
  80. Convenience overload for single axis.
  81. @param lower lowest bound that should be kept.
  82. @param upper highest bound that should be kept. If upper is inside bin interval, the
  83. whole interval is removed.
  84. */
  85. inline auto shrink(double lower, double upper) { return shrink(0, lower, upper); }
  86. /**
  87. Convenience overload for single axis.
  88. @param merge how many adjacent bins to merge into one.
  89. */
  90. inline auto rebin(unsigned merge) { return rebin(0, merge); }
  91. /**
  92. Shrink and/or rebin axes of a histogram.
  93. Returns the reduced copy of the histogram.
  94. @param hist original histogram.
  95. @param options iterable sequence of reduce_options, generated by shrink_and_rebin(),
  96. shrink(), and rebin().
  97. */
  98. #ifndef BOOST_HISTOGRAM_DOXYGEN_INVOKED
  99. template <class Histogram, class Iterable, class = detail::requires_iterable<Iterable>>
  100. #else
  101. template <class Histogram, class Iterable>
  102. #endif
  103. decltype(auto) reduce(const Histogram& hist, const Iterable& options) {
  104. const auto& old_axes = unsafe_access::axes(hist);
  105. struct option_item : reduce_option {
  106. int begin, end;
  107. };
  108. auto opts = detail::make_stack_buffer<option_item>(old_axes);
  109. for (const auto& o : options) {
  110. auto& oi = opts[o.iaxis];
  111. if (oi.merge > 0) // did we already set the option for this axis?
  112. BOOST_THROW_EXCEPTION(std::invalid_argument("indices must be unique"));
  113. oi.lower = o.lower;
  114. oi.upper = o.upper;
  115. oi.merge = o.merge;
  116. }
  117. // make new axes container with default-constructed axis instances
  118. auto axes = detail::make_default(old_axes);
  119. detail::static_if<detail::is_tuple<decltype(axes)>>(
  120. [](auto&, const auto&) {},
  121. [](auto& axes, const auto& old_axes) {
  122. axes.reserve(old_axes.size());
  123. for (const auto& a : old_axes) axes.emplace_back(detail::make_default(a));
  124. },
  125. axes, old_axes);
  126. // override default-constructed axis instances with modified instances
  127. unsigned iaxis = 0;
  128. hist.for_each_axis([&](const auto& a) {
  129. using T = detail::remove_cvref_t<decltype(a)>;
  130. auto& o = opts[iaxis];
  131. o.begin = 0;
  132. o.end = a.size();
  133. if (o.merge > 0) { // option is set?
  134. if (o.lower < o.upper) {
  135. while (o.begin != o.end && a.value(o.begin) < o.lower) ++o.begin;
  136. while (o.end != o.begin && a.value(o.end - 1) >= o.upper) --o.end;
  137. } else if (o.lower > o.upper) {
  138. // for inverted axis::regular
  139. while (o.begin != o.end && a.value(o.begin) > o.lower) ++o.begin;
  140. while (o.end != o.begin && a.value(o.end - 1) <= o.upper) --o.end;
  141. }
  142. o.end -= (o.end - o.begin) % o.merge;
  143. auto a2 = T(a, o.begin, o.end, o.merge);
  144. axis::get<T>(detail::axis_get(axes, iaxis)) = a2;
  145. } else {
  146. o.merge = 1;
  147. axis::get<T>(detail::axis_get(axes, iaxis)) = a;
  148. }
  149. ++iaxis;
  150. });
  151. auto storage = detail::make_default(unsafe_access::storage(hist));
  152. auto result = Histogram(std::move(axes), std::move(storage));
  153. auto idx = detail::make_stack_buffer<int>(unsafe_access::axes(result));
  154. for (auto x : indexed(hist, coverage::all)) {
  155. auto i = idx.begin();
  156. auto o = opts.begin();
  157. for (auto j : x.indices()) {
  158. *i = (j - o->begin);
  159. if (*i <= -1)
  160. *i = -1;
  161. else {
  162. *i /= o->merge;
  163. const int end = (o->end - o->begin) / o->merge;
  164. if (*i > end) *i = end;
  165. }
  166. ++i;
  167. ++o;
  168. }
  169. result.at(idx) += *x;
  170. }
  171. return result;
  172. }
  173. /**
  174. Shrink and/or rebin axes of a histogram.
  175. Returns the modified copy.
  176. @param hist original histogram.
  177. @param opt reduce option generated by shrink_and_rebin(), shrink(), and rebin().
  178. @param opts more reduce_options.
  179. */
  180. template <class Histogram, class... Ts>
  181. decltype(auto) reduce(const Histogram& hist, const reduce_option& opt, Ts&&... opts) {
  182. // this must be in one line, because any of the ts could be a temporary
  183. return reduce(hist, std::initializer_list<reduce_option>{opt, opts...});
  184. }
  185. } // namespace algorithm
  186. } // namespace histogram
  187. } // namespace boost
  188. #endif