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/assert.hpp>
  11 #include <boost/histogram/axis/traits.hpp>
  12 #include <boost/histogram/detail/axes.hpp>
  13 #include <boost/histogram/detail/make_default.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 <cmath>
  20 #include <initializer_list>
  21 #include <stdexcept>
  22 #include <string>
  23
  24 namespace boost {
  25 namespace histogram {
  26 namespace detail {
  27 struct reduce_option {
  28   unsigned iaxis = 0;
  29   bool indices_set = false;
  30   axis::index_type begin = 0, end = 0;
  31   bool values_set = false;
  32   double lower = 0.0, upper = 0.0;
  33   unsigned merge = 0;
  34 };
  35 } // namespace detail
  36
  37 namespace algorithm {
  38
  39 using reduce_option = detail::reduce_option;
  40
  41 /**
  42   Shrink and rebin option to be used in reduce().
  43
  44   To shrink and rebin in one command. Equivalent to passing both the shrink() and the
  45   rebin() option for the same axis to reduce.
  46
  47   @param iaxis which axis to operate on.
  48   @param lower lowest bound that should be kept.
  49   @param upper highest bound that should be kept. If upper is inside bin interval, the
  50   whole interval is removed.
  51   @param merge how many adjacent bins to merge into one.
  52  */
  53 inline reduce_option shrink_and_rebin(unsigned iaxis, double lower, double upper,
  54                                       unsigned merge) {
  55   if (lower == upper)
  56     BOOST_THROW_EXCEPTION(std::invalid_argument("lower != upper required"));
  57   if (merge == 0) BOOST_THROW_EXCEPTION(std::invalid_argument("merge > 0 required"));
  58   return {iaxis, false, 0, 0, true, lower, upper, merge};
  59 }
  60
  61 /**
  62   Slice and rebin option to be used in reduce().
  63
  64   To slice and rebin in one command. Equivalent to passing both the slice() and the
  65   rebin() option for the same axis to reduce.
  66
  67   @param iaxis which axis to operate on.
  68   @param begin first index that should be kept.
  69   @param end one past the last index that should be kept.
  70   @param merge how many adjacent bins to merge into one.
  71  */
  72 inline reduce_option slice_and_rebin(unsigned iaxis, axis::index_type begin,
  73                                      axis::index_type end, unsigned merge) {
  74   if (!(begin < end))
  75     BOOST_THROW_EXCEPTION(std::invalid_argument("begin < end required"));
  76   if (merge == 0) BOOST_THROW_EXCEPTION(std::invalid_argument("merge > 0 required"));
  77   return {iaxis, true, begin, end, false, 0.0, 0.0, merge};
  78 }
  79
  80 /**
  81   Shrink option to be used in reduce().
  82
  83   The shrink is inclusive. The bin which contains the first value starts the range of bins
  84   to keep. The bin which contains the second value is the last included in that range.
  85   When the second value is exactly equal to a lower bin edge, then the previous bin is
  86   the last in the range.
  87
  88   @param iaxis which axis to operate on.
  89   @param lower bin which contains lower is first to be kept.
  90   @param upper bin which contains upper is last to be kept, except if upper is equal to
  91   the lower edge.
  92  */
  93 inline reduce_option shrink(unsigned iaxis, double lower, double upper) {
  94   return shrink_and_rebin(iaxis, lower, upper, 1);
  95 }
  96
  97 /**
  98   Slice option to be used in reduce().
  99
 100   @param iaxis which axis to operate on.
 101   @param begin first index that should be kept.
 102   @param end one past the last index that should be kept.
 103  */
 104 inline reduce_option slice(unsigned iaxis, axis::index_type begin, axis::index_type end) {
 105   return slice_and_rebin(iaxis, begin, end, 1);
 106 }
 107
 108 /**
 109   Rebin option to be used in reduce().
 110
 111   @param iaxis which axis to operate on.
 112   @param merge how many adjacent bins to merge into one.
 113  */
 114 inline reduce_option rebin(unsigned iaxis, unsigned merge) {
 115   if (merge == 0) BOOST_THROW_EXCEPTION(std::invalid_argument("merge > 0 required"));
 116   return reduce_option{iaxis, false, 0, 0, false, 0.0, 0.0, merge};
 117 }
 118
 119 /**
 120   Shrink and rebin option to be used in reduce() (convenience overload for
 121   single axis).
 122
 123   @param lower lowest bound that should be kept.
 124   @param upper highest bound that should be kept. If upper is inside bin interval, the
 125   whole interval is removed.
 126   @param merge how many adjacent bins to merge into one.
 127 */
 128 inline reduce_option shrink_and_rebin(double lower, double upper, unsigned merge) {
 129   return shrink_and_rebin(0, lower, upper, merge);
 130 }
 131
 132 /**
 133   Slice and rebin option to be used in reduce() (convenience for 1D histograms).
 134
 135   @param begin first index that should be kept.
 136   @param end one past the last index that should be kept.
 137   @param merge how many adjacent bins to merge into one.
 138 */
 139 inline reduce_option slice_and_rebin(axis::index_type begin, axis::index_type end,
 140                                      unsigned merge) {
 141   return slice_and_rebin(0, begin, end, merge);
 142 }
 143
 144 /**
 145   Shrink option to be used in reduce() (convenience for 1D histograms).
 146
 147   @param lower lowest bound that should be kept.
 148   @param upper highest bound that should be kept. If upper is inside bin interval, the
 149   whole interval is removed.
 150 */
 151 inline reduce_option shrink(double lower, double upper) {
 152   return shrink(0, lower, upper);
 153 }
 154
 155 /**
 156   Slice option to be used in reduce() (convenience for 1D histograms).
 157
 158   @param begin first index that should be kept.
 159   @param end one past the last index that should be kept.
 160 */
 161 inline reduce_option slice(axis::index_type begin, axis::index_type end) {
 162   return slice(0, begin, end);
 163 }
 164
 165 /**
 166   Rebin option to be used in reduce() (convenience for 1D histograms).
 167
 168   @param merge how many adjacent bins to merge into one.
 169 */
 170 inline reduce_option rebin(unsigned merge) { return rebin(0, merge); }
 171
 172 /**
 173   Shrink, slice, and/or rebin axes of a histogram.
 174
 175   Returns the reduced copy of the histogram.
 176
 177   Shrinking only works with axes that accept double values. Some axis types do not support
 178   the reduce operation, for example, the builtin category axis, which is not ordered.
 179   Custom axis types must implement a special constructor (see concepts) to be reducible.
 180
 181   @param hist original histogram.
 182   @param options iterable sequence of reduce options, generated by shrink_and_rebin(),
 183   slice_and_rebin(), shrink(), slice(), and rebin().
 184  */
 185 template <class Histogram, class Iterable, class = detail::requires_iterable<Iterable>>
 186 decltype(auto) reduce(const Histogram& hist, const Iterable& options) {
 187   const auto& old_axes = unsafe_access::axes(hist);
 188
 189   auto opts = detail::make_stack_buffer<reduce_option>(old_axes);
 190   for (const reduce_option& o_in : options) {
 191     BOOST_ASSERT(o_in.merge > 0);
 192     if (o_in.iaxis >= hist.rank())
 193       BOOST_THROW_EXCEPTION(std::invalid_argument("invalid axis index"));
 194     reduce_option& o_out = opts[o_in.iaxis];
 195     if (o_out.merge > 0) {
 196       // some option was already set for this axis, see if we can merge requests
 197       if (o_in.merge > 1 && o_out.merge > 1)
 198         BOOST_THROW_EXCEPTION(std::invalid_argument("conflicting merge requests"));
 199       if ((o_in.indices_set || o_in.values_set) &&
 200           (o_out.indices_set || o_out.values_set))
 201         BOOST_THROW_EXCEPTION(
 202             std::invalid_argument("conflicting slice or shrink requests"));
 203     }
 204     if (o_in.values_set) {
 205       o_out.values_set = true;
 206       o_out.lower = o_in.lower;
 207       o_out.upper = o_in.upper;
 208     } else if (o_in.indices_set) {
 209       o_out.indices_set = true;
 210       o_out.begin = o_in.begin;
 211       o_out.end = o_in.end;
 212     }
 213     o_out.merge = (std::max)(o_in.merge, o_out.merge);
 214   }
 215
 216   // make new axes container with default-constructed axis instances
 217   auto axes = detail::make_default(old_axes);
 218   detail::static_if<detail::is_tuple<decltype(axes)>>(
 219       [](auto&, const auto&) {},
 220       [](auto& axes, const auto& old_axes) {
 221         axes.reserve(old_axes.size());
 222         detail::for_each_axis(old_axes, [&axes](const auto& a) {
 223           axes.emplace_back(detail::make_default(a));
 224         });
 225       },
 226       axes, old_axes);
 227
 228   // override default-constructed axis instances with modified instances
 229   unsigned iaxis = 0;
 230   hist.for_each_axis([&](const auto& a) {
 231     using A = std::decay_t<decltype(a)>;
 232     auto& o = opts[iaxis];
 233     if (o.merge > 0) { // option is set?
 234       detail::static_if_c<axis::traits::is_reducible<A>::value>(
 235           [&o](auto&& aout, const auto& ain) {
 236             using A = std::decay_t<decltype(ain)>;
 237             if (!o.indices_set && !o.values_set) {
 238               o.begin = 0;
 239               o.end = ain.size();
 240             } else {
 241               if (o.values_set) {
 242                 o.begin = axis::traits::index(ain, o.lower);
 243                 o.end = axis::traits::index(ain, o.upper);
 244                 if (axis::traits::value_as<double>(ain, o.end) != o.upper) ++o.end;
 245               }
 246               o.begin = (std::max)(0, o.begin);
 247               o.end = (std::min)(o.end, ain.size());
 248             }
 249             o.end -= (o.end - o.begin) % o.merge;
 250             aout = A(ain, o.begin, o.end, o.merge);
 251           },
 252           [iaxis](auto&&, const auto&) {
 253             BOOST_THROW_EXCEPTION(std::invalid_argument("axis " + std::to_string(iaxis) +
 254                                                         " is not reducible"));
 255           },
 256           axis::get<A>(detail::axis_get(axes, iaxis)), a);
 257     } else {
 258       o.merge = 1;
 259       o.begin = 0;
 260       o.end = a.size();
 261       axis::get<A>(detail::axis_get(axes, iaxis)) = a;
 262     }
 263     ++iaxis;
 264   });
 265
 266   auto storage = detail::make_default(unsafe_access::storage(hist));
 267   auto result = Histogram(std::move(axes), std::move(storage));
 268
 269   auto idx = detail::make_stack_buffer<int>(unsafe_access::axes(result));
 270   for (auto&& x : indexed(hist, coverage::all)) {
 271     auto i = idx.begin();
 272     auto o = opts.begin();
 273     for (auto j : x.indices()) {
 274       *i = (j - o->begin);
 275       if (*i <= -1)
 276         *i = -1;
 277       else {
 278         *i /= o->merge;
 279         const int end = (o->end - o->begin) / o->merge;
 280         if (*i > end) *i = end;
 281       }
 282       ++i;
 283       ++o;
 284     }
 285     result.at(idx) += *x;
 286   }
 287
 288   return result;
 289 }
 290
 291 /**
 292   Shrink, slice, and/or rebin axes of a histogram.
 293
 294   Returns the reduced copy of the histogram.
 295
 296   Shrinking only works with axes that accept double values. Some axis types do not support
 297   the reduce operation, for example, the builtin category axis, which is not ordered.
 298   Custom axis types must implement a special constructor (see concepts) to be reducible.
 299
 300   @param hist original histogram.
 301   @param opt  reduce option generated by shrink_and_rebin(), shrink(), and rebin().
 302   @param opts more reduce options.
 303  */
 304 template <class Histogram, class... Ts>
 305 decltype(auto) reduce(const Histogram& hist, const reduce_option& opt,
 306                       const Ts&... opts) {
 307   // this must be in one line, because any of the ts could be a temporary
 308   return reduce(hist, std::initializer_list<reduce_option>{opt, opts...});
 309 }
 310
 311 } // namespace algorithm
 312 } // namespace histogram
 313 } // namespace boost
 314
 315 #endif