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 */
shrink(unsigned iaxis,double lower,double upper)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 */
shrink(double lower,double upper)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 */
crop(unsigned iaxis,double lower,double upper)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 */
crop(double lower,double upper)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 */
slice(unsigned iaxis,axis::index_type begin,axis::index_type end,slice_mode mode=slice_mode::shrink)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 */
slice(axis::index_type begin,axis::index_type end,slice_mode mode=slice_mode::shrink)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 */
rebin(unsigned iaxis,unsigned merge)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 */
rebin(unsigned merge)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 */
shrink_and_rebin(unsigned iaxis,double lower,double upper,unsigned merge)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 */
shrink_and_rebin(double lower,double upper,unsigned merge)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 */
crop_and_rebin(unsigned iaxis,double lower,double upper,unsigned merge)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 */
crop_and_rebin(double lower,double upper,unsigned merge)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 */
slice_and_rebin(unsigned iaxis,axis::index_type begin,axis::index_type end,unsigned merge,slice_mode mode=slice_mode::shrink)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 */
slice_and_rebin(axis::index_type begin,axis::index_type end,unsigned merge,slice_mode mode=slice_mode::shrink)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 =
__anon5de85a100102(std::size_t iaxis, const auto& a_in) 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>
reduce(const Histogram & hist,const reduce_command & opt,const Ts &...opts)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
454