histogram/algorithm/reduce.hpp

0001 // Copyright 2018-2019 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_REDUCE_HPP
0008 #define BOOST_HISTOGRAM_ALGORITHM_REDUCE_HPP
0009
0010 #include <boost/histogram/axis/traits.hpp>
0011 #include <boost/histogram/detail/axes.hpp>
0012 #include <boost/histogram/detail/make_default.hpp>
0013 #include <boost/histogram/detail/reduce_command.hpp>
0014 #include <boost/histogram/detail/static_if.hpp>
0015 #include <boost/histogram/fwd.hpp>
0016 #include <boost/histogram/indexed.hpp>
0017 #include <boost/histogram/unsafe_access.hpp>
0018 #include <boost/throw_exception.hpp>
0019 #include <cassert>
0020 #include <cmath>
0021 #include <initializer_list>
0022 #include <stdexcept>
0023 #include <string>
0024
0025 namespace boost {
0026 namespace histogram {
0027 namespace algorithm {
0028
0029 /** Holder for a reduce command.
0030
0031   Use this type to store reduce commands in a container. The internals of this type are an
0032   implementation detail.
0033 */
0034 using reduce_command = detail::reduce_command;
0035
0036 /** Shrink command to be used in `reduce`.
0037
0038   Command is applied to axis with given index.
0039
0040   Shrinking is based on an inclusive value interval. The bin which contains the first
0041   value starts the range of bins to keep. The bin which contains the second value is the
0042   last included in that range. When the second value is exactly equal to a lower bin edge,
0043   then the previous bin is the last in the range.
0044
0045   The counts in removed bins are added to the corresponding underflow and overflow bins,
0046   if they are present. If they are not present, the counts are discarded. Also see
0047   `crop`, which always discards the counts.
0048
0049   @param iaxis which axis to operate on.
0050   @param lower bin which contains lower is first to be kept.
0051   @param upper bin which contains upper is last to be kept, except if upper is equal to
0052   the lower edge.
0053 */
0054 inline reduce_command shrink(unsigned iaxis, double lower, double upper) {
0055   if (lower == upper)
0056     BOOST_THROW_EXCEPTION(std::invalid_argument("lower != upper required"));
0057   reduce_command r;
0058   r.iaxis = iaxis;
0059   r.range = reduce_command::range_t::values;
0060   r.begin.value = lower;
0061   r.end.value = upper;
0062   r.merge = 1;
0063   r.crop = false;
0064   return r;
0065 }
0066
0067 /** Shrink command to be used in `reduce`.
0068
0069   Command is applied to corresponding axis in order of reduce arguments.
0070
0071   Shrinking is based on an inclusive value interval. The bin which contains the first
0072   value starts the range of bins to keep. The bin which contains the second value is the
0073   last included in that range. When the second value is exactly equal to a lower bin edge,
0074   then the previous bin is the last in the range.
0075
0076   The counts in removed bins are added to the corresponding underflow and overflow bins,
0077   if they are present. If they are not present, the counts are discarded. Also see
0078   `crop`, which always discards the counts.
0079
0080   @param lower bin which contains lower is first to be kept.
0081   @param upper bin which contains upper is last to be kept, except if upper is equal to
0082   the lower edge.
0083 */
0084 inline reduce_command shrink(double lower, double upper) {
0085   return shrink(reduce_command::unset, lower, upper);
0086 }
0087
0088 /** Crop command to be used in `reduce`.
0089
0090   Command is applied to axis with given index.
0091
0092   Works like `shrink` (see shrink documentation for details), but counts in removed
0093   bins are always discarded, whether underflow and overflow bins are present or not.
0094
0095   @param iaxis which axis to operate on.
0096   @param lower bin which contains lower is first to be kept.
0097   @param upper bin which contains upper is last to be kept, except if upper is equal to
0098   the lower edge.
0099 */
0100 inline reduce_command crop(unsigned iaxis, double lower, double upper) {
0101   reduce_command r = shrink(iaxis, lower, upper);
0102   r.crop = true;
0103   return r;
0104 }
0105
0106 /** Crop command to be used in `reduce`.
0107
0108   Command is applied to corresponding axis in order of reduce arguments.
0109
0110   Works like `shrink` (see shrink documentation for details), but counts in removed bins
0111   are discarded, whether underflow and overflow bins are present or not. If the cropped
0112   range goes beyond the axis range, then the content of the underflow
0113   or overflow bin which overlaps with the range is kept.
0114
0115   If the counts in an existing underflow or overflow bin are discared by the crop, the
0116   corresponding memory cells are not physically removed. Only their contents are set to
0117   zero. This technical limitation may be lifted in the future, then crop may completely
0118   remove the cropped memory cells.
0119
0120   @param lower bin which contains lower is first to be kept.
0121   @param upper bin which contains upper is last to be kept, except if upper is equal to
0122   the lower edge.
0123 */
0124 inline reduce_command crop(double lower, double upper) {
0125   return crop(reduce_command::unset, lower, upper);
0126 }
0127
0128 ///  Whether to behave like `shrink` or `crop` regarding removed bins.
0129 enum class slice_mode { shrink, crop };
0130
0131 /** Slice command to be used in `reduce`.
0132
0133   Command is applied to axis with given index.
0134
0135   Slicing works like `shrink` or `crop`, but uses bin indices instead of values.
0136
0137   @param iaxis which axis to operate on.
0138   @param begin first index that should be kept.
0139   @param end one past the last index that should be kept.
0140   @param mode whether to behave like `shrink` or `crop` regarding removed bins.
0141 */
0142 inline reduce_command slice(unsigned iaxis, axis::index_type begin, axis::index_type end,
0143                             slice_mode mode = slice_mode::shrink) {
0144   if (!(begin < end))
0145     BOOST_THROW_EXCEPTION(std::invalid_argument("begin < end required"));
0146
0147   reduce_command r;
0148   r.iaxis = iaxis;
0149   r.range = reduce_command::range_t::indices;
0150   r.begin.index = begin;
0151   r.end.index = end;
0152   r.merge = 1;
0153   r.crop = mode == slice_mode::crop;
0154   return r;
0155 }
0156
0157 /** Slice command to be used in `reduce`.
0158
0159   Command is applied to corresponding axis in order of reduce arguments.
0160
0161   Slicing works like `shrink` or `crop`, but uses bin indices instead of values.
0162
0163   @param begin first index that should be kept.
0164   @param end one past the last index that should be kept.
0165   @param mode whether to behave like `shrink` or `crop` regarding removed bins.
0166 */
0167 inline reduce_command slice(axis::index_type begin, axis::index_type end,
0168                             slice_mode mode = slice_mode::shrink) {
0169   return slice(reduce_command::unset, begin, end, mode);
0170 }
0171
0172 /** Rebin command to be used in `reduce`.
0173
0174   Command is applied to axis with given index.
0175
0176   The command merges N adjacent bins into one. This makes the axis coarser and the bins
0177   wider. The original number of bins is divided by N. If there is a rest to this devision,
0178   the axis is implicitly shrunk at the upper end by that rest.
0179
0180   @param iaxis which axis to operate on.
0181   @param merge how many adjacent bins to merge into one.
0182 */
0183 inline reduce_command rebin(unsigned iaxis, unsigned merge) {
0184   if (merge == 0) BOOST_THROW_EXCEPTION(std::invalid_argument("merge > 0 required"));
0185   reduce_command r;
0186   r.iaxis = iaxis;
0187   r.merge = merge;
0188   r.range = reduce_command::range_t::none;
0189   r.crop = false;
0190   return r;
0191 }
0192
0193 /** Rebin command to be used in `reduce`.
0194
0195   Command is applied to corresponding axis in order of reduce arguments.
0196
0197   The command merges N adjacent bins into one. This makes the axis coarser and the bins
0198   wider. The original number of bins is divided by N. If there is a rest to this devision,
0199   the axis is implicitly shrunk at the upper end by that rest.
0200
0201   @param merge how many adjacent bins to merge into one.
0202 */
0203 inline reduce_command rebin(unsigned merge) {
0204   return rebin(reduce_command::unset, merge);
0205 }
0206
0207 /** Shrink and rebin command to be used in `reduce`.
0208
0209   Command is applied to corresponding axis in order of reduce arguments.
0210
0211   To shrink(unsigned, double, double) and rebin(unsigned, unsigned) in one command (see
0212   the respective commands for more details). Equivalent to passing both commands for the
0213   same axis to `reduce`.
0214
0215   @param iaxis which axis to operate on.
0216   @param lower lowest bound that should be kept.
0217   @param upper highest bound that should be kept. If upper is inside bin interval, the
0218   whole interval is removed.
0219   @param merge how many adjacent bins to merge into one.
0220 */
0221 inline reduce_command shrink_and_rebin(unsigned iaxis, double lower, double upper,
0222                                        unsigned merge) {
0223   reduce_command r = shrink(iaxis, lower, upper);
0224   r.merge = rebin(merge).merge;
0225   return r;
0226 }
0227
0228 /** Shrink and rebin command to be used in `reduce`.
0229
0230   Command is applied to corresponding axis in order of reduce arguments.
0231
0232   To `shrink` and `rebin` in one command (see the respective commands for more
0233   details). Equivalent to passing both commands for the same axis to `reduce`.
0234
0235   @param lower lowest bound that should be kept.
0236   @param upper highest bound that should be kept. If upper is inside bin interval, the
0237   whole interval is removed.
0238   @param merge how many adjacent bins to merge into one.
0239 */
0240 inline reduce_command shrink_and_rebin(double lower, double upper, unsigned merge) {
0241   return shrink_and_rebin(reduce_command::unset, lower, upper, merge);
0242 }
0243
0244 /** Crop and rebin command to be used in `reduce`.
0245
0246   Command is applied to axis with given index.
0247
0248   To `crop` and `rebin` in one command (see the respective commands for more
0249   details). Equivalent to passing both commands for the same axis to `reduce`.
0250
0251   @param iaxis which axis to operate on.
0252   @param lower lowest bound that should be kept.
0253   @param upper highest bound that should be kept. If upper is inside bin interval,
0254   the whole interval is removed.
0255   @param merge how many adjacent bins to merge into one.
0256 */
0257 inline reduce_command crop_and_rebin(unsigned iaxis, double lower, double upper,
0258                                      unsigned merge) {
0259   reduce_command r = crop(iaxis, lower, upper);
0260   r.merge = rebin(merge).merge;
0261   return r;
0262 }
0263
0264 /** Crop and rebin command to be used in `reduce`.
0265
0266   Command is applied to corresponding axis in order of reduce arguments.
0267
0268   To `crop` and `rebin` in one command (see the respective commands for more
0269   details). Equivalent to passing both commands for the same axis to `reduce`.
0270
0271   @param lower lowest bound that should be kept.
0272   @param upper highest bound that should be kept. If upper is inside bin interval,
0273   the whole interval is removed.
0274   @param merge how many adjacent bins to merge into one.
0275 */
0276 inline reduce_command crop_and_rebin(double lower, double upper, unsigned merge) {
0277   return crop_and_rebin(reduce_command::unset, lower, upper, merge);
0278 }
0279
0280 /** Slice and rebin command to be used in `reduce`.
0281
0282   Command is applied to axis with given index.
0283
0284   To `slice` and `rebin` in one command (see the respective commands for more
0285   details). Equivalent to passing both commands for the same axis to `reduce`.
0286
0287   @param iaxis which axis to operate on.
0288   @param begin first index that should be kept.
0289   @param end one past the last index that should be kept.
0290   @param merge how many adjacent bins to merge into one.
0291   @param mode slice mode, see slice_mode.
0292 */
0293 inline reduce_command slice_and_rebin(unsigned iaxis, axis::index_type begin,
0294                                       axis::index_type end, unsigned merge,
0295                                       slice_mode mode = slice_mode::shrink) {
0296   reduce_command r = slice(iaxis, begin, end, mode);
0297   r.merge = rebin(merge).merge;
0298   return r;
0299 }
0300
0301 /** Slice and rebin command to be used in `reduce`.
0302
0303   Command is applied to corresponding axis in order of reduce arguments.
0304
0305   To `slice` and `rebin` in one command (see the respective commands for more
0306   details). Equivalent to passing both commands for the same axis to `reduce`.
0307
0308   @param begin first index that should be kept.
0309   @param end one past the last index that should be kept.
0310   @param merge how many adjacent bins to merge into one.
0311   @param mode slice mode, see slice_mode.
0312 */
0313 inline reduce_command slice_and_rebin(axis::index_type begin, axis::index_type end,
0314                                       unsigned merge,
0315                                       slice_mode mode = slice_mode::shrink) {
0316   return slice_and_rebin(reduce_command::unset, begin, end, merge, mode);
0317 }
0318
0319 /** Shrink, crop, slice, and/or rebin axes of a histogram.
0320
0321   Returns a new reduced histogram and leaves the original histogram untouched.
0322
0323   The commands `rebin` and `shrink` or `slice` for the same axis are
0324   automatically combined, this is not an error. Passing a `shrink` and a `slice`
0325   command for the same axis or two `rebin` commands triggers an `invalid_argument`
0326   exception. Trying to reducing a non-reducible axis triggers an `invalid_argument`
0327   exception. Histograms with  non-reducible axes can still be reduced along the
0328   other axes that are reducible.
0329
0330   An overload allows one to pass reduce_command as positional arguments.
0331
0332   @param hist original histogram.
0333   @param options iterable sequence of reduce commands: `shrink`, `slice`, `rebin`,
0334   `shrink_and_rebin`, or `slice_and_rebin`. The element type of the iterable should be
0335   `reduce_command`.
0336 */
0337 template <class Histogram, class Iterable, class = detail::requires_iterable<Iterable>>
0338 Histogram reduce(const Histogram& hist, const Iterable& options) {
0339   using axis::index_type;
0340
0341   const auto& old_axes = unsafe_access::axes(hist);
0342   auto opts = detail::make_stack_buffer(old_axes, reduce_command{});
0343   detail::normalize_reduce_commands(opts, options);
0344
0345   auto axes =
0346       detail::axes_transform(old_axes, [&opts](std::size_t iaxis, const auto& a_in) {
0347         using A = std::decay_t<decltype(a_in)>;
0348         using AO = axis::traits::get_options<A>;
0349         auto& o = opts[iaxis];
0350         o.is_ordered = axis::traits::ordered(a_in);
0351         if (o.merge > 0) { // option is set?
0352           o.use_underflow_bin = AO::test(axis::option::underflow);
0353           o.use_overflow_bin = AO::test(axis::option::overflow);
0354           return detail::static_if_c<axis::traits::is_reducible<A>::value>(
0355               [&o](const auto& a_in) {
0356                 if (o.range == reduce_command::range_t::none) {
0357                   // no range restriction, pure rebin
0358                   o.begin.index = 0;
0359                   o.end.index = a_in.size();
0360                 } else {
0361                   // range striction, convert values to indices as needed
0362                   if (o.range == reduce_command::range_t::values) {
0363                     const auto end_value = o.end.value;
0364                     o.begin.index = axis::traits::index(a_in, o.begin.value);
0365                     o.end.index = axis::traits::index(a_in, o.end.value);
0366                     // end = index + 1, unless end_value equal to upper bin edge
0367                     if (axis::traits::value_as<double>(a_in, o.end.index) != end_value)
0368                       ++o.end.index;
0369                   }
0370
0371                   // crop flow bins if index range does not include them
0372                   if (o.crop) {
0373                     o.use_underflow_bin &= o.begin.index < 0;
0374                     o.use_overflow_bin &= o.end.index > a_in.size();
0375                   }
0376
0377                   // now limit [begin, end] to [0, size()]
0378                   if (o.begin.index < 0) o.begin.index = 0;
0379                   if (o.end.index > a_in.size()) o.end.index = a_in.size();
0380                 }
0381                 // shorten the index range to a multiple of o.merge;
0382                 // example [1, 4] with merge = 2 is reduced to [1, 3]
0383                 o.end.index -=
0384                     (o.end.index - o.begin.index) % static_cast<index_type>(o.merge);
0385                 using A = std::decay_t<decltype(a_in)>;
0386                 return A(a_in, o.begin.index, o.end.index, o.merge);
0387               },
0388               [iaxis](const auto& a_in) {
0389                 return BOOST_THROW_EXCEPTION(std::invalid_argument(
0390                            "axis " + std::to_string(iaxis) + " is not reducible")),
0391                        a_in;
0392               },
0393               a_in);
0394         } else {
0395           // command was not set for this axis; fill noop values and copy original axis
0396           o.use_underflow_bin = AO::test(axis::option::underflow);
0397           o.use_overflow_bin = AO::test(axis::option::overflow);
0398           o.merge = 1;
0399           o.begin.index = 0;
0400           o.end.index = a_in.size();
0401           return a_in;
0402         }
0403       });
0404
0405   auto result =
0406       Histogram(std::move(axes), detail::make_default(unsafe_access::storage(hist)));
0407
0408   auto idx = detail::make_stack_buffer<index_type>(unsafe_access::axes(result));
0409   for (auto&& x : indexed(hist, coverage::all)) {
0410     auto i = idx.begin();
0411     auto o = opts.begin();
0412     bool skip = false;
0413
0414     for (auto j : x.indices()) {
0415       *i = (j - o->begin.index);
0416       if (o->is_ordered && *i <= -1) {
0417         *i = -1;
0418         if (!o->use_underflow_bin) skip = true;
0419       } else {
0420         if (*i >= 0)
0421           *i /= static_cast<index_type>(o->merge);
0422         else
0423           *i = o->end.index;
0424         const auto reduced_axis_end =
0425             (o->end.index - o->begin.index) / static_cast<index_type>(o->merge);
0426         if (*i >= reduced_axis_end) {
0427           *i = reduced_axis_end;
0428           if (!o->use_overflow_bin) skip = true;
0429         }
0430       }
0431
0432       ++i;
0433       ++o;
0434     }
0435
0436     if (!skip) result.at(idx) += *x;
0437   }
0438
0439   return result;
0440 }
0441
0442 /** Shrink, slice, and/or rebin axes of a histogram.
0443
0444   Returns a new reduced histogram and leaves the original histogram untouched.
0445
0446   The commands `rebin` and `shrink` or `slice` for the same axis are
0447   automatically combined, this is not an error. Passing a `shrink` and a `slice`
0448   command for the same axis or two `rebin` commands triggers an invalid_argument
0449   exception. It is safe to reduce histograms with some axis that are not reducible along
0450   the other axes. Trying to reducing a non-reducible axis triggers an invalid_argument
0451   exception.
0452
0453   An overload allows one to pass an iterable of reduce_command.
0454
0455   @param hist original histogram.
0456   @param opt first reduce command; one of `shrink`, `slice`, `rebin`,
0457   `shrink_and_rebin`, or `slice_or_rebin`.
0458   @param opts more reduce commands.
0459 */
0460 template <class Histogram, class... Ts>
0461 Histogram reduce(const Histogram& hist, const reduce_command& opt, const Ts&... opts) {
0462   // this must be in one line, because any of the ts could be a temporary
0463   return reduce(hist, std::initializer_list<reduce_command>{opt, opts...});
0464 }
0465
0466 } // namespace algorithm
0467 } // namespace histogram
0468 } // namespace boost
0469
0470 #endif