Skip to content

Commit bbf4104

Browse files
hahnjohahnjo
committed
[hist] Implement RHistEngine::Slice
It puts together RSliceSpec, axis slicing, and RSliceBinIndexMapper to transform the bin contents and return the sliced histogram.
1 parent cc83bbc commit bbf4104

2 files changed

Lines changed: 337 additions & 0 deletions

File tree

‎hist/histv7/inc/ROOT/RHistEngine.hxx‎

Lines changed: 130 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -12,6 +12,8 @@
1212
#include "RHistUtils.hxx"
1313
#include "RLinearizedIndex.hxx"
1414
#include "RRegularAxis.hxx"
15+
#include "RSliceBinIndexMapper.hxx"
16+
#include "RSliceSpec.hxx"
1517
#include "RWeight.hxx"
1618

1719
#include <array>
@@ -639,6 +641,134 @@ public:
639641
}
640642
}
641643

644+
/// Slice this histogram with an RSliceSpec per dimension.
645+
///
646+
/// With a range, only the specified bins are retained. All other bin contents are transferred to the underflow and
647+
/// overflow bins:
648+
/// \code
649+
/// ROOT::Experimental::RHistEngine<int> hist(/* one dimension */);
650+
/// // Fill the histogram with a number of entries...
651+
/// auto sliced = hist.Slice({hist.GetAxes()[0].GetNormalRange(1, 5)});
652+
/// // The returned histogram will have 4 normal bins, an underflow and an overflow bin.
653+
/// \endcode
654+
///
655+
/// Slicing can also perform operations per dimension, see RSliceSpec. RSliceSpec::ROperationRebin allows to rebin
656+
/// the histogram axis, grouping a number of normal bins into a new one:
657+
/// \code
658+
/// ROOT::Experimental::RHistEngine<int> hist(/* one dimension */);
659+
/// // Fill the histogram with a number of entries...
660+
/// auto rebinned = hist.Slice(ROOT::Experimental::RSliceSpec::ROperationRebin(2));
661+
/// // The returned histogram has groups of two normal bins merged.
662+
/// \endcode
663+
///
664+
/// RSliceSpec::ROperationSum sums the bin contents along that axis, which allows to project to a lower-dimensional
665+
/// histogram:
666+
/// \code
667+
/// ROOT::Experimental::RHistEngine<int> hist({/* two dimensions */});
668+
/// // Fill the histogram with a number of entries...
669+
/// auto projected = hist.Slice(ROOT::Experimental::RSliceSpec{}, ROOT::Experimental::RSliceSpec::ROperationSum{});
670+
/// // The returned histogram has one dimension, with bin contents summed along the second axis.
671+
/// \endcode
672+
/// Note that it is not allowed to sum along all histogram axes because the return value would be a scalar.
673+
///
674+
/// Ranges and operations can be combined. In that case, the range is applied before the operation.
675+
///
676+
/// \param[in] sliceSpecs the slice specifications for each axis
677+
/// \return the sliced histogram
678+
/// \par See also
679+
/// the \ref Slice(const A &... args) const "variadic function template overload" accepting arguments directly
680+
RHistEngine Slice(const std::vector<RSliceSpec> &sliceSpecs) const
681+
{
682+
if (sliceSpecs.size() != GetNDimensions()) {
683+
throw std::invalid_argument("invalid number of specifications passed to Slice");
684+
}
685+
686+
// Slice the axes.
687+
std::vector<RAxisVariant> axes;
688+
for (std::size_t i = 0; i < sliceSpecs.size(); i++) {
689+
// A sum operation makes the dimension disappear.
690+
if (sliceSpecs[i].GetOperationSum() == nullptr) {
691+
axes.push_back(fAxes.Get()[i].Slice(sliceSpecs[i]));
692+
}
693+
}
694+
if (axes.empty()) {
695+
throw std::invalid_argument("summing across all dimensions is not supported");
696+
}
697+
698+
RHistEngine sliced(std::move(axes));
699+
700+
// Create the helper objects to map the bin contents to the sliced histogram.
701+
Internal::RSliceBinIndexMapper mapper(sliceSpecs);
702+
assert(mapper.GetMappedDimensionality() == sliced.GetNDimensions());
703+
std::vector<RBinIndex> mappedIndices(mapper.GetMappedDimensionality());
704+
705+
auto origRange = fAxes.GetFullMultiDimRange();
706+
auto origRangeIt = origRange.begin();
707+
708+
for (std::size_t i = 0; i < fBinContents.size(); i++) {
709+
const auto &origIndices = *origRangeIt;
710+
#ifndef NDEBUG
711+
// Verify that the original indices correspond to the iteration variable.
712+
RLinearizedIndex origIndex = fAxes.ComputeGlobalIndex(origIndices);
713+
assert(origIndex.fValid);
714+
assert(origIndex.fIndex == i);
715+
#endif
716+
717+
bool success = mapper.Map(origIndices, mappedIndices);
718+
if (success) {
719+
RLinearizedIndex mappedIndex = sliced.fAxes.ComputeGlobalIndex(mappedIndices);
720+
assert(mappedIndex.fValid);
721+
sliced.fBinContents[mappedIndex.fIndex] += fBinContents[i];
722+
}
723+
++origRangeIt;
724+
}
725+
726+
return sliced;
727+
}
728+
729+
/// Slice this histogram with an RSliceSpec per dimension.
730+
///
731+
/// With a range, only the specified bins are retained. All other bin contents are transferred to the underflow and
732+
/// overflow bins:
733+
/// \code
734+
/// ROOT::Experimental::RHistEngine<int> hist(/* one dimension */);
735+
/// // Fill the histogram with a number of entries...
736+
/// auto sliced = hist.Slice(hist.GetAxes()[0].GetNormalRange(1, 5));
737+
/// // The returned histogram will have 4 normal bins, an underflow and an overflow bin.
738+
/// \endcode
739+
///
740+
/// Slicing can also perform operations per dimension, see RSliceSpec. RSliceSpec::ROperationRebin allows to rebin
741+
/// the histogram axis, grouping a number of normal bins into a new one:
742+
/// \code
743+
/// ROOT::Experimental::RHistEngine<int> hist(/* one dimension */);
744+
/// // Fill the histogram with a number of entries...
745+
/// auto rebinned = hist.Slice(ROOT::Experimental::RSliceSpec::ROperationRebin(2));
746+
/// // The returned histogram has groups of two normal bins merged.
747+
/// \endcode
748+
///
749+
/// RSliceSpec::ROperationSum sums the bin contents along that axis, which allows to project to a lower-dimensional
750+
/// histogram:
751+
/// \code
752+
/// ROOT::Experimental::RHistEngine<int> hist({/* two dimensions */});
753+
/// // Fill the histogram with a number of entries...
754+
/// auto projected = hist.Slice(ROOT::Experimental::RSliceSpec{}, ROOT::Experimental::RSliceSpec::ROperationSum{});
755+
/// // The returned histogram has one dimension, with bin contents summed along the second axis.
756+
/// \endcode
757+
/// Note that it is not allowed to sum along all histogram axes because the return value would be a scalar.
758+
///
759+
/// Ranges and operations can be combined. In that case, the range is applied before the operation.
760+
///
761+
/// \param[in] args the arguments for each axis
762+
/// \return the sliced histogram
763+
/// \par See also
764+
/// the \ref Slice(const std::vector<RSliceSpec> &sliceSpecs) const "function overload" accepting `std::vector`
765+
template <typename... A>
766+
RHistEngine Slice(const A &...args) const
767+
{
768+
std::vector<RSliceSpec> sliceSpecs{args...};
769+
return Slice(sliceSpecs);
770+
}
771+
642772
/// %ROOT Streamer function to throw when trying to store an object of this class.
643773
void Streamer(TBuffer &) { throw std::runtime_error("unable to store RHistEngine"); }
644774
};

‎hist/histv7/test/hist_slice.cxx‎

Lines changed: 207 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -285,3 +285,210 @@ TEST(RSliceBinIndexMapper, MapSliceFullSum)
285285
EXPECT_TRUE(success);
286286
}
287287
}
288+
289+
TEST(RHistEngine, SliceInvalidNumberOfArguments)
290+
{
291+
static constexpr std::size_t Bins = 20;
292+
const RRegularAxis axis(Bins, {0, Bins});
293+
const RHistEngine<int> engine1(axis);
294+
ASSERT_EQ(engine1.GetNDimensions(), 1);
295+
const RHistEngine<int> engine2(axis, axis);
296+
ASSERT_EQ(engine2.GetNDimensions(), 2);
297+
298+
EXPECT_NO_THROW(engine1.Slice(RSliceSpec{}));
299+
EXPECT_THROW(engine1.Slice(RSliceSpec{}, RSliceSpec{}), std::invalid_argument);
300+
301+
EXPECT_THROW(engine2.Slice(RSliceSpec{}), std::invalid_argument);
302+
EXPECT_NO_THROW(engine2.Slice(RSliceSpec{}, RSliceSpec{}));
303+
EXPECT_THROW(engine2.Slice(RSliceSpec{}, RSliceSpec{}, RSliceSpec{}), std::invalid_argument);
304+
}
305+
306+
TEST(RHistEngine, SliceSumAll)
307+
{
308+
static constexpr std::size_t Bins = 20;
309+
const RRegularAxis axis(Bins, {0, Bins});
310+
const RHistEngine<int> engine1(axis);
311+
ASSERT_EQ(engine1.GetNDimensions(), 1);
312+
const RHistEngine<int> engine2(axis, axis);
313+
ASSERT_EQ(engine2.GetNDimensions(), 2);
314+
315+
const RSliceSpec sum(RSliceSpec::ROperationSum{});
316+
EXPECT_THROW(engine1.Slice(sum), std::invalid_argument);
317+
EXPECT_THROW(engine2.Slice(sum, sum), std::invalid_argument);
318+
}
319+
320+
TEST(RHistEngine, SliceFull)
321+
{
322+
static constexpr std::size_t Bins = 20;
323+
const RRegularAxis axis(Bins, {0, Bins});
324+
RHistEngine<int> engine(axis);
325+
326+
engine.SetBinContent(RBinIndex::Underflow(), 100);
327+
for (std::size_t i = 0; i < Bins; i++) {
328+
engine.SetBinContent(i, i + 1);
329+
}
330+
engine.SetBinContent(RBinIndex::Overflow(), 200);
331+
332+
// Three different ways of "slicing" which will keep the entire axis.
333+
for (auto sliceSpec : {RSliceSpec{}, RSliceSpec(axis.GetFullRange()), RSliceSpec(axis.GetNormalRange())}) {
334+
const auto sliced = engine.Slice(sliceSpec);
335+
ASSERT_EQ(sliced.GetNDimensions(), 1);
336+
EXPECT_TRUE(sliced.GetAxes()[0].GetRegularAxis() != nullptr);
337+
EXPECT_EQ(sliced.GetAxes()[0].GetNNormalBins(), Bins);
338+
EXPECT_EQ(sliced.GetTotalNBins(), Bins + 2);
339+
340+
EXPECT_EQ(sliced.GetBinContent(RBinIndex::Underflow()), 100);
341+
for (std::size_t i = 0; i < Bins; i++) {
342+
EXPECT_EQ(sliced.GetBinContent(i), i + 1);
343+
}
344+
EXPECT_EQ(sliced.GetBinContent(RBinIndex::Overflow()), 200);
345+
}
346+
}
347+
348+
TEST(RHistEngine, SliceNormal)
349+
{
350+
static constexpr std::size_t Bins = 20;
351+
const RRegularAxis axis(Bins, {0, Bins});
352+
RHistEngine<int> engine(axis);
353+
354+
engine.SetBinContent(RBinIndex::Underflow(), 100);
355+
for (std::size_t i = 0; i < Bins; i++) {
356+
engine.SetBinContent(i, i + 1);
357+
}
358+
engine.SetBinContent(RBinIndex::Overflow(), 200);
359+
360+
const auto sliced = engine.Slice(axis.GetNormalRange(1, 5));
361+
ASSERT_EQ(sliced.GetNDimensions(), 1);
362+
const auto *regular = sliced.GetAxes()[0].GetRegularAxis();
363+
ASSERT_TRUE(regular != nullptr);
364+
EXPECT_EQ(regular->GetNNormalBins(), 4);
365+
EXPECT_EQ(regular->GetLow(), 1);
366+
EXPECT_EQ(regular->GetHigh(), 5);
367+
EXPECT_EQ(sliced.GetTotalNBins(), 6);
368+
369+
EXPECT_EQ(sliced.GetBinContent(RBinIndex::Underflow()), 101);
370+
for (std::size_t i = 0; i < 4; i++) {
371+
EXPECT_EQ(sliced.GetBinContent(i), i + 2);
372+
}
373+
EXPECT_EQ(sliced.GetBinContent(RBinIndex::Overflow()), 395);
374+
}
375+
376+
TEST(RHistEngine, SliceRebin)
377+
{
378+
static constexpr std::size_t Bins = 20;
379+
const RRegularAxis axis(Bins, {0, Bins});
380+
RHistEngine<int> engine(axis);
381+
382+
engine.SetBinContent(RBinIndex::Underflow(), 100);
383+
for (std::size_t i = 0; i < Bins; i++) {
384+
engine.SetBinContent(i, i + 1);
385+
}
386+
engine.SetBinContent(RBinIndex::Overflow(), 200);
387+
388+
const auto sliced = engine.Slice(RSliceSpec::ROperationRebin(2));
389+
ASSERT_EQ(sliced.GetNDimensions(), 1);
390+
const auto *regular = sliced.GetAxes()[0].GetRegularAxis();
391+
ASSERT_TRUE(regular != nullptr);
392+
EXPECT_EQ(regular->GetNNormalBins(), Bins / 2);
393+
EXPECT_EQ(regular->GetLow(), 0);
394+
EXPECT_EQ(regular->GetHigh(), Bins);
395+
EXPECT_EQ(sliced.GetTotalNBins(), Bins / 2 + 2);
396+
397+
EXPECT_EQ(sliced.GetBinContent(RBinIndex::Underflow()), 100);
398+
for (std::size_t i = 0; i < Bins / 2; i++) {
399+
EXPECT_EQ(sliced.GetBinContent(i), 4 * i + 3);
400+
}
401+
EXPECT_EQ(sliced.GetBinContent(RBinIndex::Overflow()), 200);
402+
}
403+
404+
TEST(RHistEngine, SliceRangeRebin)
405+
{
406+
static constexpr std::size_t Bins = 20;
407+
const RRegularAxis axis(Bins, {0, Bins});
408+
RHistEngine<int> engine(axis);
409+
410+
engine.SetBinContent(RBinIndex::Underflow(), 100);
411+
for (std::size_t i = 0; i < Bins; i++) {
412+
engine.SetBinContent(i, i + 1);
413+
}
414+
engine.SetBinContent(RBinIndex::Overflow(), 200);
415+
416+
const RSliceSpec spec(axis.GetNormalRange(1, 5), RSliceSpec::ROperationRebin(2));
417+
const auto sliced = engine.Slice(spec);
418+
ASSERT_EQ(sliced.GetNDimensions(), 1);
419+
const auto *regular = sliced.GetAxes()[0].GetRegularAxis();
420+
ASSERT_TRUE(regular != nullptr);
421+
EXPECT_EQ(regular->GetNNormalBins(), 2);
422+
EXPECT_EQ(regular->GetLow(), 1);
423+
EXPECT_EQ(regular->GetHigh(), 5);
424+
EXPECT_EQ(sliced.GetTotalNBins(), 4);
425+
426+
EXPECT_EQ(sliced.GetBinContent(RBinIndex::Underflow()), 101);
427+
for (std::size_t i = 0; i < 2; i++) {
428+
EXPECT_EQ(sliced.GetBinContent(i), 4 * i + 5);
429+
}
430+
EXPECT_EQ(sliced.GetBinContent(RBinIndex::Overflow()), 395);
431+
}
432+
433+
TEST(RHistEngine, SliceSum)
434+
{
435+
static constexpr std::size_t Bins = 20;
436+
const RRegularAxis axis(Bins, {0, Bins});
437+
RHistEngine<int> engine(axis, axis);
438+
439+
engine.SetBinContent(RBinIndex::Underflow(), 0, 1000);
440+
engine.SetBinContent(RBinIndex::Underflow(), 2, 2000);
441+
for (std::size_t i = 0; i < Bins; i++) {
442+
for (std::size_t j = 0; j < Bins; j++) {
443+
engine.SetBinContent(i, RBinIndex::Underflow(), 100 * i);
444+
engine.SetBinContent(i, RBinIndex(j), i * Bins + j);
445+
engine.SetBinContent(i, RBinIndex::Overflow(), 200 * i);
446+
}
447+
}
448+
engine.SetBinContent(RBinIndex::Overflow(), 3, 3000);
449+
engine.SetBinContent(RBinIndex::Overflow(), 6, 4000);
450+
451+
const auto sliced = engine.Slice(RSliceSpec{}, RSliceSpec::ROperationSum{});
452+
ASSERT_EQ(sliced.GetNDimensions(), 1);
453+
EXPECT_TRUE(sliced.GetAxes()[0].GetRegularAxis() != nullptr);
454+
EXPECT_EQ(sliced.GetAxes()[0].GetNNormalBins(), Bins);
455+
EXPECT_EQ(sliced.GetTotalNBins(), Bins + 2);
456+
457+
EXPECT_EQ(sliced.GetBinContent(RBinIndex::Underflow()), 3000);
458+
for (std::size_t i = 0; i < Bins; i++) {
459+
EXPECT_EQ(sliced.GetBinContent(i), i * (100 + Bins * Bins + 200) + Bins * (Bins - 1) / 2);
460+
}
461+
EXPECT_EQ(sliced.GetBinContent(RBinIndex::Overflow()), 7000);
462+
}
463+
464+
TEST(RHistEngine, SliceRangeSum)
465+
{
466+
static constexpr std::size_t Bins = 20;
467+
const RRegularAxis axis(Bins, {0, Bins});
468+
RHistEngine<int> engine(axis, axis);
469+
470+
engine.SetBinContent(RBinIndex::Underflow(), 0, 1000);
471+
engine.SetBinContent(RBinIndex::Underflow(), 2, 2000);
472+
for (std::size_t i = 0; i < Bins; i++) {
473+
for (std::size_t j = 0; j < Bins; j++) {
474+
engine.SetBinContent(i, RBinIndex::Underflow(), 100 * i);
475+
engine.SetBinContent(i, RBinIndex(j), i * Bins + j);
476+
engine.SetBinContent(i, RBinIndex::Overflow(), 200 * i);
477+
}
478+
}
479+
engine.SetBinContent(RBinIndex::Overflow(), 3, 3000);
480+
engine.SetBinContent(RBinIndex::Overflow(), 6, 4000);
481+
482+
const RSliceSpec rangeSum(axis.GetNormalRange(1, 5), RSliceSpec::ROperationSum{});
483+
const auto sliced = engine.Slice(RSliceSpec{}, rangeSum);
484+
ASSERT_EQ(sliced.GetNDimensions(), 1);
485+
EXPECT_TRUE(sliced.GetAxes()[0].GetRegularAxis() != nullptr);
486+
EXPECT_EQ(sliced.GetAxes()[0].GetNNormalBins(), Bins);
487+
EXPECT_EQ(sliced.GetTotalNBins(), Bins + 2);
488+
489+
EXPECT_EQ(sliced.GetBinContent(RBinIndex::Underflow()), 2000);
490+
for (std::size_t i = 0; i < Bins; i++) {
491+
EXPECT_EQ(sliced.GetBinContent(i), 4 * i * Bins + 10);
492+
}
493+
EXPECT_EQ(sliced.GetBinContent(RBinIndex::Overflow()), 3000);
494+
}

0 commit comments

Comments
 (0)