| 123456789101112131415161718192021222324252627282930313233343536373839404142434445464748495051525354555657585960616263646566676869707172737475767778798081828384858687888990919293949596979899100101102103104105106107108109110111112113114115116117118119120121122123124125126127128129130131132133134135136137138139140141142143144145146147148149150151152153154155156157158159160161162163164165166167168169170171172173174175176177178179180181182183184185186187188189190191192193194195196197198199200201202203204205206207208209210211212213214215216217218219220221222223 |
- // Copyright 2018 Hans Dembinski
- //
- // Distributed under the Boost Software License, Version 1.0.
- // (See accompanying file LICENSE_1_0.txt
- // or copy at http://www.boost.org/LICENSE_1_0.txt)
- #ifndef BOOST_HISTOGRAM_ALGORITHM_REDUCE_HPP
- #define BOOST_HISTOGRAM_ALGORITHM_REDUCE_HPP
- #include <boost/assert.hpp>
- #include <boost/histogram/detail/axes.hpp>
- #include <boost/histogram/detail/meta.hpp>
- #include <boost/histogram/fwd.hpp>
- #include <boost/histogram/indexed.hpp>
- #include <boost/histogram/unsafe_access.hpp>
- #include <boost/throw_exception.hpp>
- #include <cmath>
- #include <limits>
- #include <stdexcept>
- #include <type_traits>
- namespace boost {
- namespace histogram {
- namespace algorithm {
- /**
- Option type returned by the helper functions shrink_and_rebin(), shrink(), rebin().
- */
- struct reduce_option {
- #ifndef BOOST_HISTOGRAM_DOXYGEN_INVOKED
- unsigned iaxis = 0;
- double lower, upper;
- unsigned merge = 0;
- reduce_option() noexcept = default;
- reduce_option(unsigned i, double l, double u, unsigned m)
- : iaxis(i), lower(l), upper(u), merge(m) {
- if (lower == upper)
- BOOST_THROW_EXCEPTION(std::invalid_argument("lower != upper required"));
- if (merge == 0) BOOST_THROW_EXCEPTION(std::invalid_argument("merge > 0 required"));
- }
- #endif
- };
- /**
- Shrink and rebin option.
- @param iaxis which axis to operate on.
- @param lower lowest bound that should be kept.
- @param upper highest bound that should be kept. If upper is inside bin interval, the
- whole interval is removed.
- @param merge how many adjacent bins to merge into one.
- */
- inline auto shrink_and_rebin(unsigned iaxis, double lower, double upper, unsigned merge) {
- return reduce_option{iaxis, lower, upper, merge};
- }
- /**
- Shrink option.
- @param iaxis which axis to operate on.
- @param lower lowest bound that should be kept.
- @param upper highest bound that should be kept. If upper is inside bin interval, the
- whole interval is removed.
- */
- inline auto shrink(unsigned iaxis, double lower, double upper) {
- return reduce_option{iaxis, lower, upper, 1};
- }
- /**
- Rebin option.
- @param iaxis which axis to operate on.
- @param merge how many adjacent bins to merge into one.
- */
- inline auto rebin(unsigned iaxis, unsigned merge) {
- return reduce_option{iaxis, std::numeric_limits<double>::quiet_NaN(),
- std::numeric_limits<double>::quiet_NaN(), merge};
- }
- /**
- Convenience overload for single axis.
- @param lower lowest bound that should be kept.
- @param upper highest bound that should be kept. If upper is inside bin interval, the
- whole interval is removed.
- @param merge how many adjacent bins to merge into one.
- */
- inline auto shrink_and_rebin(double lower, double upper, unsigned merge) {
- return shrink_and_rebin(0, lower, upper, merge);
- }
- /**
- Convenience overload for single axis.
- @param lower lowest bound that should be kept.
- @param upper highest bound that should be kept. If upper is inside bin interval, the
- whole interval is removed.
- */
- inline auto shrink(double lower, double upper) { return shrink(0, lower, upper); }
- /**
- Convenience overload for single axis.
- @param merge how many adjacent bins to merge into one.
- */
- inline auto rebin(unsigned merge) { return rebin(0, merge); }
- /**
- Shrink and/or rebin axes of a histogram.
- Returns the reduced copy of the histogram.
- @param hist original histogram.
- @param options iterable sequence of reduce_options, generated by shrink_and_rebin(),
- shrink(), and rebin().
- */
- #ifndef BOOST_HISTOGRAM_DOXYGEN_INVOKED
- template <class Histogram, class Iterable, class = detail::requires_iterable<Iterable>>
- #else
- template <class Histogram, class Iterable>
- #endif
- decltype(auto) reduce(const Histogram& hist, const Iterable& options) {
- const auto& old_axes = unsafe_access::axes(hist);
- struct option_item : reduce_option {
- int begin, end;
- };
- auto opts = detail::make_stack_buffer<option_item>(old_axes);
- for (const auto& o : options) {
- auto& oi = opts[o.iaxis];
- if (oi.merge > 0) // did we already set the option for this axis?
- BOOST_THROW_EXCEPTION(std::invalid_argument("indices must be unique"));
- oi.lower = o.lower;
- oi.upper = o.upper;
- oi.merge = o.merge;
- }
- // make new axes container with default-constructed axis instances
- auto axes = detail::make_default(old_axes);
- detail::static_if<detail::is_tuple<decltype(axes)>>(
- [](auto&, const auto&) {},
- [](auto& axes, const auto& old_axes) {
- axes.reserve(old_axes.size());
- for (const auto& a : old_axes) axes.emplace_back(detail::make_default(a));
- },
- axes, old_axes);
- // override default-constructed axis instances with modified instances
- unsigned iaxis = 0;
- hist.for_each_axis([&](const auto& a) {
- using T = detail::remove_cvref_t<decltype(a)>;
- auto& o = opts[iaxis];
- o.begin = 0;
- o.end = a.size();
- if (o.merge > 0) { // option is set?
- if (o.lower < o.upper) {
- while (o.begin != o.end && a.value(o.begin) < o.lower) ++o.begin;
- while (o.end != o.begin && a.value(o.end - 1) >= o.upper) --o.end;
- } else if (o.lower > o.upper) {
- // for inverted axis::regular
- while (o.begin != o.end && a.value(o.begin) > o.lower) ++o.begin;
- while (o.end != o.begin && a.value(o.end - 1) <= o.upper) --o.end;
- }
- o.end -= (o.end - o.begin) % o.merge;
- auto a2 = T(a, o.begin, o.end, o.merge);
- axis::get<T>(detail::axis_get(axes, iaxis)) = a2;
- } else {
- o.merge = 1;
- axis::get<T>(detail::axis_get(axes, iaxis)) = a;
- }
- ++iaxis;
- });
- auto storage = detail::make_default(unsafe_access::storage(hist));
- auto result = Histogram(std::move(axes), std::move(storage));
- auto idx = detail::make_stack_buffer<int>(unsafe_access::axes(result));
- for (auto x : indexed(hist, coverage::all)) {
- auto i = idx.begin();
- auto o = opts.begin();
- for (auto j : x.indices()) {
- *i = (j - o->begin);
- if (*i <= -1)
- *i = -1;
- else {
- *i /= o->merge;
- const int end = (o->end - o->begin) / o->merge;
- if (*i > end) *i = end;
- }
- ++i;
- ++o;
- }
- result.at(idx) += *x;
- }
- return result;
- }
- /**
- Shrink and/or rebin axes of a histogram.
- Returns the modified copy.
- @param hist original histogram.
- @param opt reduce option generated by shrink_and_rebin(), shrink(), and rebin().
- @param opts more reduce_options.
- */
- template <class Histogram, class... Ts>
- decltype(auto) reduce(const Histogram& hist, const reduce_option& opt, Ts&&... opts) {
- // this must be in one line, because any of the ts could be a temporary
- return reduce(hist, std::initializer_list<reduce_option>{opt, opts...});
- }
- } // namespace algorithm
- } // namespace histogram
- } // namespace boost
- #endif
|