ceph/src/boost/boost/histogram/algorithm/reduce.hpp

   1 // Copyright 2018-2019 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
   7 #ifndef BOOST_HISTOGRAM_ALGORITHM_REDUCE_HPP
   8 #define BOOST_HISTOGRAM_ALGORITHM_REDUCE_HPP
   9
  10 #include <boost/histogram/axis/traits.hpp>
  11 #include <boost/histogram/detail/axes.hpp>
  12 #include <boost/histogram/detail/make_default.hpp>
  13 #include <boost/histogram/detail/reduce_command.hpp>
  14 #include <boost/histogram/detail/static_if.hpp>
  15 #include <boost/histogram/fwd.hpp>
  16 #include <boost/histogram/indexed.hpp>
  17 #include <boost/histogram/unsafe_access.hpp>
  18 #include <boost/throw_exception.hpp>
  19 #include <cassert>
  20 #include <cmath>
  21 #include <initializer_list>
  22 #include <stdexcept>
  23 #include <string>
  24
  25 namespace boost {
  26 namespace histogram {
  27 namespace algorithm {
  28
  29 /** Holder for a reduce command.
  30
  31   Use this type to store reduce commands in a container. The internals of this type are an
  32   implementation detail.
  33 */
  34 using reduce_command = detail::reduce_command;
  35
  36 using reduce_option [[deprecated("use reduce_command instead")]] =
  37     reduce_command; ///< deprecated
  38
  39 /** Shrink command to be used in `reduce`.
  40
  41   Command is applied to axis with given index.
  42
  43   Shrinking is based on an inclusive value interval. The bin which contains the first
  44   value starts the range of bins to keep. The bin which contains the second value is the
  45   last included in that range. When the second value is exactly equal to a lower bin edge,
  46   then the previous bin is the last in the range.
  47
  48   The counts in removed bins are added to the corresponding underflow and overflow bins,
  49   if they are present. If they are not present, the counts are discarded. Also see
  50   `crop`, which always discards the counts.
  51
  52   @param iaxis which axis to operate on.
  53   @param lower bin which contains lower is first to be kept.
  54   @param upper bin which contains upper is last to be kept, except if upper is equal to
  55   the lower edge.
  56 */
  57 inline reduce_command shrink(unsigned iaxis, double lower, double upper) {
  58   if (lower == upper)
  59     BOOST_THROW_EXCEPTION(std::invalid_argument("lower != upper required"));
  60   reduce_command r;
  61   r.iaxis = iaxis;
  62   r.range = reduce_command::range_t::values;
  63   r.begin.value = lower;
  64   r.end.value = upper;
  65   r.merge = 1;
  66   r.crop = false;
  67   return r;
  68 }
  69
  70 /** Shrink command to be used in `reduce`.
  71
  72   Command is applied to corresponding axis in order of reduce arguments.
  73
  74   Shrinking is based on an inclusive value interval. The bin which contains the first
  75   value starts the range of bins to keep. The bin which contains the second value is the
  76   last included in that range. When the second value is exactly equal to a lower bin edge,
  77   then the previous bin is the last in the range.
  78
  79   The counts in removed bins are added to the corresponding underflow and overflow bins,
  80   if they are present. If they are not present, the counts are discarded. Also see
  81   `crop`, which always discards the counts.
  82
  83   @param lower bin which contains lower is first to be kept.
  84   @param upper bin which contains upper is last to be kept, except if upper is equal to
  85   the lower edge.
  86 */
  87 inline reduce_command shrink(double lower, double upper) {
  88   return shrink(reduce_command::unset, lower, upper);
  89 }
  90
  91 /** Crop command to be used in `reduce`.
  92
  93   Command is applied to axis with given index.
  94
  95   Works like `shrink` (see shrink documentation for details), but counts in removed
  96   bins are always discarded, whether underflow and overflow bins are present or not.
  97
  98   @param iaxis which axis to operate on.
  99   @param lower bin which contains lower is first to be kept.
 100   @param upper bin which contains upper is last to be kept, except if upper is equal to
 101   the lower edge.
 102 */
 103 inline reduce_command crop(unsigned iaxis, double lower, double upper) {
 104   reduce_command r = shrink(iaxis, lower, upper);
 105   r.crop = true;
 106   return r;
 107 }
 108
 109 /** Crop command to be used in `reduce`.
 110
 111   Command is applied to corresponding axis in order of reduce arguments.
 112
 113   Works like `shrink` (see shrink documentation for details), but counts in removed bins
 114   are discarded, whether underflow and overflow bins are present or not.
 115
 116   @param lower bin which contains lower is first to be kept.
 117   @param upper bin which contains upper is last to be kept, except if upper is equal to
 118   the lower edge.
 119 */
 120 inline reduce_command crop(double lower, double upper) {
 121   return crop(reduce_command::unset, lower, upper);
 122 }
 123
 124 ///  Whether to behave like `shrink` or `crop` regarding removed bins.
 125 enum class slice_mode { shrink, crop };
 126
 127 /** Slice command to be used in `reduce`.
 128
 129   Command is applied to axis with given index.
 130
 131   Slicing works like `shrink` or `crop`, but uses bin indices instead of values.
 132
 133   @param iaxis which axis to operate on.
 134   @param begin first index that should be kept.
 135   @param end one past the last index that should be kept.
 136   @param mode whether to behave like `shrink` or `crop` regarding removed bins.
 137 */
 138 inline reduce_command slice(unsigned iaxis, axis::index_type begin, axis::index_type end,
 139                             slice_mode mode = slice_mode::shrink) {
 140   if (!(begin < end))
 141     BOOST_THROW_EXCEPTION(std::invalid_argument("begin < end required"));
 142
 143   reduce_command r;
 144   r.iaxis = iaxis;
 145   r.range = reduce_command::range_t::indices;
 146   r.begin.index = begin;
 147   r.end.index = end;
 148   r.merge = 1;
 149   r.crop = mode == slice_mode::crop;
 150   return r;
 151 }
 152
 153 /** Slice command to be used in `reduce`.
 154
 155   Command is applied to corresponding axis in order of reduce arguments.
 156
 157   Slicing works like `shrink` or `crop`, but uses bin indices instead of values.
 158
 159   @param begin first index that should be kept.
 160   @param end one past the last index that should be kept.
 161   @param mode whether to behave like `shrink` or `crop` regarding removed bins.
 162 */
 163 inline reduce_command slice(axis::index_type begin, axis::index_type end,
 164                             slice_mode mode = slice_mode::shrink) {
 165   return slice(reduce_command::unset, begin, end, mode);
 166 }
 167
 168 /** Rebin command to be used in `reduce`.
 169
 170   Command is applied to axis with given index.
 171
 172   The command merges N adjacent bins into one. This makes the axis coarser and the bins
 173   wider. The original number of bins is divided by N. If there is a rest to this devision,
 174   the axis is implicitly shrunk at the upper end by that rest.
 175
 176   @param iaxis which axis to operate on.
 177   @param merge how many adjacent bins to merge into one.
 178 */
 179 inline reduce_command rebin(unsigned iaxis, unsigned merge) {
 180   if (merge == 0) BOOST_THROW_EXCEPTION(std::invalid_argument("merge > 0 required"));
 181   reduce_command r;
 182   r.iaxis = iaxis;
 183   r.merge = merge;
 184   r.range = reduce_command::range_t::none;
 185   r.crop = false;
 186   return r;
 187 }
 188
 189 /** Rebin command to be used in `reduce`.
 190
 191   Command is applied to corresponding axis in order of reduce arguments.
 192
 193   The command merges N adjacent bins into one. This makes the axis coarser and the bins
 194   wider. The original number of bins is divided by N. If there is a rest to this devision,
 195   the axis is implicitly shrunk at the upper end by that rest.
 196
 197   @param merge how many adjacent bins to merge into one.
 198 */
 199 inline reduce_command rebin(unsigned merge) {
 200   return rebin(reduce_command::unset, merge);
 201 }
 202
 203 /** Shrink and rebin command to be used in `reduce`.
 204
 205   Command is applied to corresponding axis in order of reduce arguments.
 206
 207   To shrink(unsigned, double, double) and rebin(unsigned, unsigned) in one command (see
 208   the respective commands for more details). Equivalent to passing both commands for the
 209   same axis to `reduce`.
 210
 211   @param iaxis which axis to operate on.
 212   @param lower lowest bound that should be kept.
 213   @param upper highest bound that should be kept. If upper is inside bin interval, the
 214   whole interval is removed.
 215   @param merge how many adjacent bins to merge into one.
 216 */
 217 inline reduce_command shrink_and_rebin(unsigned iaxis, double lower, double upper,
 218                                        unsigned merge) {
 219   reduce_command r = shrink(iaxis, lower, upper);
 220   r.merge = rebin(merge).merge;
 221   return r;
 222 }
 223
 224 /** Shrink and rebin command to be used in `reduce`.
 225
 226   Command is applied to corresponding axis in order of reduce arguments.
 227
 228   To `shrink` and `rebin` in one command (see the respective commands for more
 229   details). Equivalent to passing both commands for the same axis to `reduce`.
 230
 231   @param lower lowest bound that should be kept.
 232   @param upper highest bound that should be kept. If upper is inside bin interval, the
 233   whole interval is removed.
 234   @param merge how many adjacent bins to merge into one.
 235 */
 236 inline reduce_command shrink_and_rebin(double lower, double upper, unsigned merge) {
 237   return shrink_and_rebin(reduce_command::unset, lower, upper, merge);
 238 }
 239
 240 /** Crop and rebin command to be used in `reduce`.
 241
 242   Command is applied to axis with given index.
 243
 244   To `crop` and `rebin` in one command (see the respective commands for more
 245   details). Equivalent to passing both commands for the same axis to `reduce`.
 246
 247   @param iaxis which axis to operate on.
 248   @param lower lowest bound that should be kept.
 249   @param upper highest bound that should be kept. If upper is inside bin interval,
 250   the whole interval is removed.
 251   @param merge how many adjacent bins to merge into one.
 252 */
 253 inline reduce_command crop_and_rebin(unsigned iaxis, double lower, double upper,
 254                                      unsigned merge) {
 255   reduce_command r = crop(iaxis, lower, upper);
 256   r.merge = rebin(merge).merge;
 257   return r;
 258 }
 259
 260 /** Crop and rebin command to be used in `reduce`.
 261
 262   Command is applied to corresponding axis in order of reduce arguments.
 263
 264   To `crop` and `rebin` in one command (see the respective commands for more
 265   details). Equivalent to passing both commands for the same axis to `reduce`.
 266
 267   @param lower lowest bound that should be kept.
 268   @param upper highest bound that should be kept. If upper is inside bin interval,
 269   the whole interval is removed.
 270   @param merge how many adjacent bins to merge into one.
 271 */
 272 inline reduce_command crop_and_rebin(double lower, double upper, unsigned merge) {
 273   return crop_and_rebin(reduce_command::unset, lower, upper, merge);
 274 }
 275
 276 /** Slice and rebin command to be used in `reduce`.
 277
 278   Command is applied to axis with given index.
 279
 280   To `slice` and `rebin` in one command (see the respective commands for more
 281   details). Equivalent to passing both commands for the same axis to `reduce`.
 282
 283   @param iaxis which axis to operate on.
 284   @param begin first index that should be kept.
 285   @param end one past the last index that should be kept.
 286   @param merge how many adjacent bins to merge into one.
 287   @param mode slice mode, see slice_mode.
 288 */
 289 inline reduce_command slice_and_rebin(unsigned iaxis, axis::index_type begin,
 290                                       axis::index_type end, unsigned merge,
 291                                       slice_mode mode = slice_mode::shrink) {
 292   reduce_command r = slice(iaxis, begin, end, mode);
 293   r.merge = rebin(merge).merge;
 294   return r;
 295 }
 296
 297 /** Slice and rebin command to be used in `reduce`.
 298
 299   Command is applied to corresponding axis in order of reduce arguments.
 300
 301   To `slice` and `rebin` in one command (see the respective commands for more
 302   details). Equivalent to passing both commands for the same axis to `reduce`.
 303
 304   @param begin first index that should be kept.
 305   @param end one past the last index that should be kept.
 306   @param merge how many adjacent bins to merge into one.
 307   @param mode slice mode, see slice_mode.
 308 */
 309 inline reduce_command slice_and_rebin(axis::index_type begin, axis::index_type end,
 310                                       unsigned merge,
 311                                       slice_mode mode = slice_mode::shrink) {
 312   return slice_and_rebin(reduce_command::unset, begin, end, merge, mode);
 313 }
 314
 315 /** Shrink, crop, slice, and/or rebin axes of a histogram.
 316
 317   Returns a new reduced histogram and leaves the original histogram untouched.
 318
 319   The commands `rebin` and `shrink` or `slice` for the same axis are
 320   automatically combined, this is not an error. Passing a `shrink` and a `slice`
 321   command for the same axis or two `rebin` commands triggers an `invalid_argument`
 322   exception. Trying to reducing a non-reducible axis triggers an `invalid_argument`
 323   exception. Histograms with  non-reducible axes can still be reduced along the
 324   other axes that are reducible.
 325
 326   @param hist original histogram.
 327   @param options iterable sequence of reduce commands: `shrink`, `slice`, `rebin`,
 328   `shrink_and_rebin`, or `slice_and_rebin`. The element type of the iterable should be
 329   `reduce_command`.
 330 */
 331 template <class Histogram, class Iterable, class = detail::requires_iterable<Iterable>>
 332 Histogram reduce(const Histogram& hist, const Iterable& options) {
 333   using axis::index_type;
 334
 335   const auto& old_axes = unsafe_access::axes(hist);
 336   auto opts = detail::make_stack_buffer(old_axes, reduce_command{});
 337   detail::normalize_reduce_commands(opts, options);
 338
 339   auto axes =
 340       detail::axes_transform(old_axes, [&opts](std::size_t iaxis, const auto& a_in) {
 341         using A = std::decay_t<decltype(a_in)>;
 342         using AO = axis::traits::get_options<A>;
 343         auto& o = opts[iaxis];
 344         o.is_ordered = axis::traits::ordered(a_in);
 345         if (o.merge > 0) { // option is set?
 346           o.use_underflow_bin = !o.crop && AO::test(axis::option::underflow);
 347           o.use_overflow_bin = !o.crop && AO::test(axis::option::overflow);
 348           return detail::static_if_c<axis::traits::is_reducible<A>::value>(
 349               [&o](const auto& a_in) {
 350                 if (o.range == reduce_command::range_t::none) {
 351                   o.begin.index = 0;
 352                   o.end.index = a_in.size();
 353                 } else {
 354                   if (o.range == reduce_command::range_t::values) {
 355                     const auto end_value = o.end.value;
 356                     o.begin.index = axis::traits::index(a_in, o.begin.value);
 357                     o.end.index = axis::traits::index(a_in, o.end.value);
 358                     // end = index + 1, unless end_value equal to upper bin edge
 359                     if (axis::traits::value_as<double>(a_in, o.end.index) != end_value)
 360                       ++o.end.index;
 361                   }
 362                   // limit [begin, end] to [0, size()]
 363                   if (o.begin.index < 0) o.begin.index = 0;
 364                   if (o.end.index > a_in.size()) o.end.index = a_in.size();
 365                 }
 366                 // shorten the index range to a multiple of o.merge;
 367                 // example [1, 4] with merge = 2 is reduced to [1, 3]
 368                 o.end.index -=
 369                     (o.end.index - o.begin.index) % static_cast<index_type>(o.merge);
 370                 using A = std::decay_t<decltype(a_in)>;
 371                 return A(a_in, o.begin.index, o.end.index, o.merge);
 372               },
 373               [iaxis](const auto& a_in) {
 374                 return BOOST_THROW_EXCEPTION(std::invalid_argument(
 375                            "axis " + std::to_string(iaxis) + " is not reducible")),
 376                        a_in;
 377               },
 378               a_in);
 379         } else {
 380           // command was not set for this axis; fill noop values and copy original axis
 381           o.use_underflow_bin = AO::test(axis::option::underflow);
 382           o.use_overflow_bin = AO::test(axis::option::overflow);
 383           o.merge = 1;
 384           o.begin.index = 0;
 385           o.end.index = a_in.size();
 386           return a_in;
 387         }
 388       });
 389
 390   auto result =
 391       Histogram(std::move(axes), detail::make_default(unsafe_access::storage(hist)));
 392
 393   auto idx = detail::make_stack_buffer<index_type>(unsafe_access::axes(result));
 394   for (auto&& x : indexed(hist, coverage::all)) {
 395     auto i = idx.begin();
 396     auto o = opts.begin();
 397     bool skip = false;
 398
 399     for (auto j : x.indices()) {
 400       *i = (j - o->begin.index);
 401       if (o->is_ordered && *i <= -1) {
 402         *i = -1;
 403         if (!o->use_underflow_bin) skip = true;
 404       } else {
 405         if (*i >= 0)
 406           *i /= static_cast<index_type>(o->merge);
 407         else
 408           *i = o->end.index;
 409         const auto reduced_axis_end =
 410             (o->end.index - o->begin.index) / static_cast<index_type>(o->merge);
 411         if (*i >= reduced_axis_end) {
 412           *i = reduced_axis_end;
 413           if (!o->use_overflow_bin) skip = true;
 414         }
 415       }
 416
 417       ++i;
 418       ++o;
 419     }
 420
 421     if (!skip) result.at(idx) += *x;
 422   }
 423
 424   return result;
 425 }
 426
 427 /** Shrink, slice, and/or rebin axes of a histogram.
 428
 429   Returns a new reduced histogram and leaves the original histogram untouched.
 430
 431   The commands `rebin` and `shrink` or `slice` for the same axis are
 432   automatically combined, this is not an error. Passing a `shrink` and a `slice`
 433   command for the same axis or two `rebin` commands triggers an invalid_argument
 434   exception. It is safe to reduce histograms with some axis that are not reducible along
 435   the other axes. Trying to reducing a non-reducible axis triggers an invalid_argument
 436   exception.
 437
 438   @param hist original histogram.
 439   @param opt first reduce command; one of `shrink`, `slice`, `rebin`,
 440   `shrink_and_rebin`, or `slice_or_rebin`.
 441   @param opts more reduce commands.
 442 */
 443 template <class Histogram, class... Ts>
 444 Histogram reduce(const Histogram& hist, const reduce_command& opt, const Ts&... opts) {
 445   // this must be in one line, because any of the ts could be a temporary
 446   return reduce(hist, std::initializer_list<reduce_command>{opt, opts...});
 447 }
 448
 449 } // namespace algorithm
 450 } // namespace histogram
 451 } // namespace boost
 452
 453 #endif