123456789101112131415161718192021222324252627282930313233343536373839404142434445464748495051525354555657585960616263646566676869707172737475767778798081828384858687888990919293949596979899100101102103104105106107108109110111112113114115116117118119120121122123124125126127128129130131132133134135136137138139140141142143144145146147148149150151152153154155156157158159160161162163164165166167168169170171172173174175176177178179180181182183184185186187188189190191192193194195196197198199200201202203204205206207208209210211212213214215216217218219220221222223224225226227228229230231232233234235236237238239240241242243244245246247248249250251252253254255256257258259260261262263264265266267268269270271272273274275276277278279280281282283284285286287288289290291292293294295296297298299300301302303304305306307308309310311312313314315 |
- // Copyright 2018-2019 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/axis/traits.hpp>
- #include <boost/histogram/detail/axes.hpp>
- #include <boost/histogram/detail/make_default.hpp>
- #include <boost/histogram/detail/static_if.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 <initializer_list>
- #include <stdexcept>
- #include <string>
- namespace boost {
- namespace histogram {
- namespace detail {
- struct reduce_option {
- unsigned iaxis = 0;
- bool indices_set = false;
- axis::index_type begin = 0, end = 0;
- bool values_set = false;
- double lower = 0.0, upper = 0.0;
- unsigned merge = 0;
- };
- } // namespace detail
- namespace algorithm {
- using reduce_option = detail::reduce_option;
- /**
- Shrink and rebin option to be used in reduce().
- To shrink and rebin in one command. Equivalent to passing both the shrink() and the
- rebin() option for the same axis to reduce.
- @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 reduce_option shrink_and_rebin(unsigned iaxis, double lower, double upper,
- unsigned merge) {
- if (lower == upper)
- BOOST_THROW_EXCEPTION(std::invalid_argument("lower != upper required"));
- if (merge == 0) BOOST_THROW_EXCEPTION(std::invalid_argument("merge > 0 required"));
- return {iaxis, false, 0, 0, true, lower, upper, merge};
- }
- /**
- Slice and rebin option to be used in reduce().
- To slice and rebin in one command. Equivalent to passing both the slice() and the
- rebin() option for the same axis to reduce.
- @param iaxis which axis to operate on.
- @param begin first index that should be kept.
- @param end one past the last index that should be kept.
- @param merge how many adjacent bins to merge into one.
- */
- inline reduce_option slice_and_rebin(unsigned iaxis, axis::index_type begin,
- axis::index_type end, unsigned merge) {
- if (!(begin < end))
- BOOST_THROW_EXCEPTION(std::invalid_argument("begin < end required"));
- if (merge == 0) BOOST_THROW_EXCEPTION(std::invalid_argument("merge > 0 required"));
- return {iaxis, true, begin, end, false, 0.0, 0.0, merge};
- }
- /**
- Shrink option to be used in reduce().
- The shrink is inclusive. The bin which contains the first value starts the range of bins
- to keep. The bin which contains the second value is the last included in that range.
- When the second value is exactly equal to a lower bin edge, then the previous bin is
- the last in the range.
- @param iaxis which axis to operate on.
- @param lower bin which contains lower is first to be kept.
- @param upper bin which contains upper is last to be kept, except if upper is equal to
- the lower edge.
- */
- inline reduce_option shrink(unsigned iaxis, double lower, double upper) {
- return shrink_and_rebin(iaxis, lower, upper, 1);
- }
- /**
- Slice option to be used in reduce().
- @param iaxis which axis to operate on.
- @param begin first index that should be kept.
- @param end one past the last index that should be kept.
- */
- inline reduce_option slice(unsigned iaxis, axis::index_type begin, axis::index_type end) {
- return slice_and_rebin(iaxis, begin, end, 1);
- }
- /**
- Rebin option to be used in reduce().
- @param iaxis which axis to operate on.
- @param merge how many adjacent bins to merge into one.
- */
- inline reduce_option rebin(unsigned iaxis, unsigned merge) {
- if (merge == 0) BOOST_THROW_EXCEPTION(std::invalid_argument("merge > 0 required"));
- return reduce_option{iaxis, false, 0, 0, false, 0.0, 0.0, merge};
- }
- /**
- Shrink and rebin option to be used in reduce() (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 reduce_option shrink_and_rebin(double lower, double upper, unsigned merge) {
- return shrink_and_rebin(0, lower, upper, merge);
- }
- /**
- Slice and rebin option to be used in reduce() (convenience for 1D histograms).
- @param begin first index that should be kept.
- @param end one past the last index that should be kept.
- @param merge how many adjacent bins to merge into one.
- */
- inline reduce_option slice_and_rebin(axis::index_type begin, axis::index_type end,
- unsigned merge) {
- return slice_and_rebin(0, begin, end, merge);
- }
- /**
- Shrink option to be used in reduce() (convenience for 1D histograms).
- @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 reduce_option shrink(double lower, double upper) {
- return shrink(0, lower, upper);
- }
- /**
- Slice option to be used in reduce() (convenience for 1D histograms).
- @param begin first index that should be kept.
- @param end one past the last index that should be kept.
- */
- inline reduce_option slice(axis::index_type begin, axis::index_type end) {
- return slice(0, begin, end);
- }
- /**
- Rebin option to be used in reduce() (convenience for 1D histograms).
- @param merge how many adjacent bins to merge into one.
- */
- inline reduce_option rebin(unsigned merge) { return rebin(0, merge); }
- /**
- Shrink, slice, and/or rebin axes of a histogram.
- Returns the reduced copy of the histogram.
- Shrinking only works with axes that accept double values. Some axis types do not support
- the reduce operation, for example, the builtin category axis, which is not ordered.
- Custom axis types must implement a special constructor (see concepts) to be reducible.
- @param hist original histogram.
- @param options iterable sequence of reduce options, generated by shrink_and_rebin(),
- slice_and_rebin(), shrink(), slice(), and rebin().
- */
- template <class Histogram, class Iterable, class = detail::requires_iterable<Iterable>>
- decltype(auto) reduce(const Histogram& hist, const Iterable& options) {
- const auto& old_axes = unsafe_access::axes(hist);
- auto opts = detail::make_stack_buffer<reduce_option>(old_axes);
- for (const reduce_option& o_in : options) {
- BOOST_ASSERT(o_in.merge > 0);
- if (o_in.iaxis >= hist.rank())
- BOOST_THROW_EXCEPTION(std::invalid_argument("invalid axis index"));
- reduce_option& o_out = opts[o_in.iaxis];
- if (o_out.merge > 0) {
- // some option was already set for this axis, see if we can merge requests
- if (o_in.merge > 1 && o_out.merge > 1)
- BOOST_THROW_EXCEPTION(std::invalid_argument("conflicting merge requests"));
- if ((o_in.indices_set || o_in.values_set) &&
- (o_out.indices_set || o_out.values_set))
- BOOST_THROW_EXCEPTION(
- std::invalid_argument("conflicting slice or shrink requests"));
- }
- if (o_in.values_set) {
- o_out.values_set = true;
- o_out.lower = o_in.lower;
- o_out.upper = o_in.upper;
- } else if (o_in.indices_set) {
- o_out.indices_set = true;
- o_out.begin = o_in.begin;
- o_out.end = o_in.end;
- }
- o_out.merge = (std::max)(o_in.merge, o_out.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());
- detail::for_each_axis(old_axes, [&axes](const auto& a) {
- 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 A = std::decay_t<decltype(a)>;
- auto& o = opts[iaxis];
- if (o.merge > 0) { // option is set?
- detail::static_if_c<axis::traits::is_reducible<A>::value>(
- [&o](auto&& aout, const auto& ain) {
- using A = std::decay_t<decltype(ain)>;
- if (!o.indices_set && !o.values_set) {
- o.begin = 0;
- o.end = ain.size();
- } else {
- if (o.values_set) {
- o.begin = axis::traits::index(ain, o.lower);
- o.end = axis::traits::index(ain, o.upper);
- if (axis::traits::value_as<double>(ain, o.end) != o.upper) ++o.end;
- }
- o.begin = (std::max)(0, o.begin);
- o.end = (std::min)(o.end, ain.size());
- }
- o.end -= (o.end - o.begin) % o.merge;
- aout = A(ain, o.begin, o.end, o.merge);
- },
- [iaxis](auto&&, const auto&) {
- BOOST_THROW_EXCEPTION(std::invalid_argument("axis " + std::to_string(iaxis) +
- " is not reducible"));
- },
- axis::get<A>(detail::axis_get(axes, iaxis)), a);
- } else {
- o.merge = 1;
- o.begin = 0;
- o.end = a.size();
- axis::get<A>(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, slice, and/or rebin axes of a histogram.
- Returns the reduced copy of the histogram.
- Shrinking only works with axes that accept double values. Some axis types do not support
- the reduce operation, for example, the builtin category axis, which is not ordered.
- Custom axis types must implement a special constructor (see concepts) to be reducible.
- @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,
- const 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
|