Skip to content

[histv7] Division for RHist. #5497

New issue

Have a question about this project? Sign up for a free GitHub account to open an issue and contact its maintainers and the community.

Sign up for GitHub

By clicking “Sign up for GitHub”, you agree to our terms of service and privacy statement. We’ll occasionally send you account related emails.

Already on GitHub? Sign in to your account

Closed
claireguyot wants to merge 6 commits into from
Closed

[histv7] Division for RHist. #5497

Show file tree
Hide file tree
Changes from all commits
Commits
Show all changes
6 commits
Select commit Hold shift + click to select a range
baed459
[histv7] NFC: removing unused enum in RHistImpl, and adjusting commen…
claireguyot Apr 29, 2020
de88bc5
[histv7] A histogram can be divided by an other. The error isn't prop…
claireguyot Apr 29, 2020
19a8766
[histv7] A histogram can be multiplied by another.
claireguyot Apr 29, 2020
a0a0817
[histv7] A separate function tests the compatibility of two axes binn…
claireguyot May 5, 2020
bfb05ab
[histv7] NFC: Comments for Divide and DivideBinomial are more accurate.
claireguyot May 5, 2020
8662c9f
[histv7] for loops in MultiplyByScalar are simpler, and result is cor…
claireguyot May 5, 2020
File filter

Filter by extension

Filter by extension
Conversations
Failed to load comments.
Loading
Jump to
Jump to file
Failed to load files.
Loading
Diff view
Diff view
100 changes: 91 additions & 9 deletions hist/histv7/inc/ROOT/RHist.hxx
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -171,6 +171,13 @@ public:
std::swap(fFillFunc, other.fFillFunc);
}

/// Multiply all stats of a histogram by a certain scalar.
template <typename T>
RHist &operator*=(T scalar) {
fImpl->GetStat().MultiplyByScalar(scalar);
return *this;
}

private:
std::unique_ptr<ImplBase_t> fImpl; ///< The actual histogram implementation.
FillFunc_t fFillFunc = nullptr; ///< Pointer to RHistImpl::Fill() member function.
Expand Down Expand Up @@ -301,15 +308,29 @@ using RH3I = RHist<3, int, RHistStatContent>;
using RH3LL = RHist<3, int64_t, RHistStatContent>;
///\}

template <int DIMENSIONS, class PRECISION,
template <int D_, class P_> class... STAT_FIRST,
template <int D_, class P_> class... STAT_SECOND>
void SameAxesOrThrow(RHist<DIMENSIONS, PRECISION, STAT_FIRST...> &first, const RHist<DIMENSIONS, PRECISION, STAT_SECOND...> &second)
{
auto& firstImpl = *first.GetImpl();
const auto& secondImpl = *second.GetImpl();
for (int dim = 0; dim < DIMENSIONS; ++dim) {
if (!firstImpl.GetAxis(dim).HasSameBinningAs(secondImpl.GetAxis(dim))) {
throw std::runtime_error("Attempted to merge two RHists with incompatible axis binning!");
}
}
}

/// Add two histograms.
///
/// This operation may currently only be performed if the two histograms have
/// the same axis configuration, use the same precision, and if `from` records
/// at least the same statistics as `to` (recording more stats is fine).
///
/// Adding histograms with incompatible axis binning will be reported at runtime
/// with an `std::runtime_error`. Insufficient statistics in the source
/// histogram will be detected at compile-time and result in a compiler error.
/// with an `std::runtime_error`. Incompatible statistics' size in the source
Copy link
Contributor

@HadrienG2 HadrienG2 May 6, 2020

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

It's not just about size, the actual constraint (as of the current code) is that "from" must record at least the same statistics as "to" or else the toImpl.GetStat().Add(fromImpl.GetStat()) call will fail because the statistics which only exist in "to" will not find their counterpart in "from" and therefore the call to the corresponding virtual statistics class method in the implementation will fail.

/// histograms will be detected at compile-time and result in a compiler error.
///
/// In the future, we may either adopt a more relaxed definition of histogram
/// addition or provide a mechanism to convert from one histogram type to
Expand All @@ -320,19 +341,80 @@ template <int DIMENSIONS, class PRECISION,
void Add(RHist<DIMENSIONS, PRECISION, STAT_TO...> &to, const RHist<DIMENSIONS, PRECISION, STAT_FROM...> &from)
{
// Enforce "same axis configuration" policy.
auto& toImpl = *to.GetImpl();
const auto& fromImpl = *from.GetImpl();
for (int dim = 0; dim < DIMENSIONS; ++dim) {
if (!toImpl.GetAxis(dim).HasSameBinningAs(fromImpl.GetAxis(dim))) {
throw std::runtime_error("Attempted to add RHists with incompatible axis binning");
}
}
SameAxesOrThrow(to, from);

// Now that we know that the two axes have the same binning, we can just add
// the statistics directly.
auto& toImpl = *to.GetImpl();
const auto& fromImpl = *from.GetImpl();
toImpl.GetStat().Add(fromImpl.GetStat());
}

/// Divide the `hist` by the `dividor`.
///
/// `hist = hist/dividor`
/// This operation may only be performed if the two histograms have
/// the same axis configuration, use the same precision, and if `dividor` records
/// at least the same statistics as `hist` (recording more stats is fine).
///
/// Dividing histograms with incompatible axis binning will be reported at runtime
/// with an `std::runtime_error`. Incompatible statistics' size in the source
/// histograms will be detected at compile-time and result in a compiler error.
///
/// In the future, we may either adopt a more relaxed definition of histogram
/// division or provide a mechanism to convert from one histogram type to
/// another. We currently favor the latter path.
///
/// NOTE : The resulting errors are calculated assuming uncorrelated histograms,
/// using Gaussian error propagation.
/// See ROOT::Experimental::DivideBinomial for binomial error propagation.
template <int DIMENSIONS, class PRECISION,
template <int D_, class P_> class... STAT_TO,
template <int D_, class P_> class... STAT_FROM>
void Divide(RHist<DIMENSIONS, PRECISION, STAT_TO...> &hist, const RHist<DIMENSIONS, PRECISION, STAT_FROM...> &dividor)
{
// Enforce "same axis configuration" policy.
SameAxesOrThrow(hist, dividor);

// Now that we know that the two axes have the same binning, we can call
// `Divide` directly on the statistics.
auto& histImpl = *hist.GetImpl();
const auto& dividingImpl = *dividor.GetImpl();
histImpl.GetStat().Divide(dividingImpl.GetStat());
}

/// Divide the `hist` by the `dividor`.
///
/// `hist = hist/dividor`
/// This operation may only be performed if the two histograms have
/// the same axis configuration, use the same precision, and if `dividor` records
/// at least the same statistics as `hist` (recording more stats is fine).
///
/// Dividing histograms with incompatible axis binning will be reported at runtime
/// with an `std::runtime_error`. Incompatible statistics' size in the source
/// histograms will be detected at compile-time and result in a compiler error.
///
/// In the future, we may either adopt a more relaxed definition of histogram
/// division or provide a mechanism to convert from one histogram type to
/// another. We currently favor the latter path.
///
/// NOTE : The resulting errors are calculated using binomial error propagation.
/// See ROOT::Experimental::Divide for Gaussian error propagation.
template <int DIMENSIONS, class PRECISION,
template <int D_, class P_> class... STAT_TO,
template <int D_, class P_> class... STAT_FROM>
void DivideBinomial(RHist<DIMENSIONS, PRECISION, STAT_TO...> &hist, const RHist<DIMENSIONS, PRECISION, STAT_FROM...> &dividor)
{
// Enforce "same axis configuration" policy.
SameAxesOrThrow(hist, dividor);

// Now that we know that the two axes have the same binning, we can call
// `Divide` directly on the statistics.
auto& histImpl = *hist.GetImpl();
const auto& dividingImpl = *dividor.GetImpl();
histImpl.GetStat().DivideBinomial(dividingImpl.GetStat());
}

} // namespace Experimental
} // namespace ROOT

Expand Down
178 changes: 173 additions & 5 deletions hist/histv7/inc/ROOT/RHistData.hxx
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -149,7 +149,7 @@ public:
/// Retrieve the under-/overflow content array (non-const).
Content_t &GetOverflowContentArray() { return fOverflowBinContent; }

/// Merge with other RHistStatContent, assuming same bin configuration.
/// Add an other RHistStatContent, assuming same bin configuration.
void Add(const RHistStatContent& other) {
assert(fBinContent.size() == other.fBinContent.size()
&& "this and other have incompatible bin configuration!");
Expand All @@ -161,6 +161,32 @@ public:
for (size_t b = 0; b < fOverflowBinContent.size(); ++b)
fOverflowBinContent[b] += other.fOverflowBinContent[b];
}

/// Multiply by a scalar.
template <typename T>
void MultiplyByScalar(const T scalar) {
for (auto &b: fBinContent)
b *= scalar;
for (auto &b: fOverflowBinContent)
b *= scalar;
}

/// Divide by an other RHistStatContent, assuming same bin configuration.
void Divide(const RHistStatContent& other) {
assert(fBinContent.size() == other.fBinContent.size());
assert(fOverflowBinContent.size() == other.fOverflowBinContent.size());
for (size_t b = 0; b < fBinContent.size(); ++b)
fBinContent[b] /= other.fBinContent[b];
for (size_t b = 0; b < fOverflowBinContent.size(); ++b)
fOverflowBinContent[b] /= other.fOverflowBinContent[b];
fEntries += other.fEntries;
}

/// Divide by an other RHistStatContent, assuming same bin configuration.
void DivideBinomial(const RHistStatContent& other) {
Divide(other);
fEntries = other.fEntries;
}
};

/**
Expand Down Expand Up @@ -201,10 +227,27 @@ public:
/// Get the sum of weights.
Weight_t GetSumOfWeights() const { return fSumWeights; }

/// Merge with other RHistStatTotalSumOfWeights data, assuming same bin configuration.
/// Add an other RHistStatTotalSumOfWeights data, assuming same bin configuration.
void Add(const RHistStatTotalSumOfWeights& other) {
fSumWeights += other.fSumWeights;
}

/// Multiply by a scalar.
template <typename T>
void MultiplyByScalar(const T scalar) {
fSumWeights *= scalar;
}

/// Divide by an other RHistStatTotalSumOfWeights data, assuming same bin configuration.
void Divide(const RHistStatTotalSumOfWeights& other) {
fSumWeights /= other.fSumWeights;
}

/// Divide by an other RHistStatTotalSumOfWeights data, assuming same bin configuration.
void DivideBinomial(const RHistStatTotalSumOfWeights& other) {
// /!\ Should be the sum of the updated weights (after division) of the bin contents
Copy link
Member

@Axel-Naumann Axel-Naumann Apr 30, 2020

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Out of curiosity: what's /!\ Looks great for a warning sign, does Doxygen understand that? And maybe:

Suggested change
// /!\ Should be the sum of the updated weights (after division) of the bin contents
// /!\ Must be called *after* the bin contents have been divided!

(if that's what it is). More occurrences below.

Divide(other);
}
};

/**
Expand Down Expand Up @@ -245,10 +288,32 @@ public:
/// Get the sum of weights.
Weight_t GetSumOfSquaredWeights() const { return fSumWeights2; }

/// Merge with other RHistStatTotalSumOfSquaredWeights data, assuming same bin configuration.
/// Add an other RHistStatTotalSumOfSquaredWeights data, assuming same bin configuration.
void Add(const RHistStatTotalSumOfSquaredWeights& other) {
fSumWeights2 += other.fSumWeights2;
}

/// Multiply by a scalar.
template <typename T>
void MultiplyByScalar(const T scalar) {
fSumWeights2 *= scalar * scalar;
}

/// Divide by an other RHistStatTotalSumOfSquaredWeights data, assuming same bin configuration.
void Divide(const RHistStatTotalSumOfSquaredWeights& other) {
fSumWeights2 /= other.fSumWeights2;
// fSumWeights2 = (fSumWeights2 * other.fUnchangedBinContent * other.fUnchangedBinContent
// + other.fSumWeights2 * fUnchangedBinContent * fUnchangedBinContent)
// / (other.fUnchangedBinContent * other.fUnchangedBinContent * other.fUnchangedBinContent * other.fUnchangedBinContent);
}

/// Divide by an other RHistStatTotalSumOfSquaredWeights data, assuming same bin configuration.
void DivideBinomial(const RHistStatTotalSumOfSquaredWeights& other) {
Divide(other);
// fSumWeights2 = std::abs(((1. - 2. * fUnchangedBinContent / other.fUnchangedBinContent) * fSumWeights2
// + fUnchangedBinContent * fUnchangedBinContent * other.fSumWeights2 / other.fUnchangedBinContent * other.fUnchangedBinContent)
// / (other.fUnchangedBinContent * other.fUnchangedBinContent));
}
};

/**
Expand Down Expand Up @@ -358,7 +423,7 @@ public:
/// Get the structure holding the under-/overflow sum of squares of weights (non-const).
std::vector<double> &GetOverflowSumOfSquaredWeights() { return fOverflowSumWeightsSquared; }

/// Merge with other `RHistStatUncertainty` data, assuming same bin configuration.
/// Add an other `RHistStatUncertainty` data, assuming same bin configuration.
void Add(const RHistStatUncertainty& other) {
assert(fSumWeightsSquared.size() == other.fSumWeightsSquared.size()
&& "this and other have incompatible bin configuration!");
Expand All @@ -369,6 +434,46 @@ public:
for (size_t b = 0; b < fOverflowSumWeightsSquared.size(); ++b)
fOverflowSumWeightsSquared[b] += other.fOverflowSumWeightsSquared[b];
}

/// Multiply by a scalar.
template <typename T>
void MultiplyByScalar(const T scalar) {
for (auto &b: fSumWeightsSquared)
b *= scalar * scalar;
for (auto &b: fOverflowSumWeightsSquared)
b *= scalar * scalar;
}

/// Divide by an other `RHistStatUncertainty` data, assuming same bin configuration.
void Divide(const RHistStatUncertainty& other) {
assert(fSumWeightsSquared.size() == other.fSumWeightsSquared.size());
assert(fOverflowSumWeightsSquared.size() == other.fOverflowSumWeightsSquared.size());
for (size_t b = 0; b < fSumWeightsSquared.size(); ++b)
fSumWeightsSquared[b] /= other.fSumWeightsSquared[b];
for (size_t b = 0; b < fOverflowSumWeightsSquared.size(); ++b)
fOverflowSumWeightsSquared[b] /= other.fOverflowSumWeightsSquared[b];
// for (size_t b = 0; b < fSumWeightsSquared.size(); ++b)
// fSumWeightsSquared[b] = (fSumWeightsSquared[b] * other.fUnchangedBinContent * other.fUnchangedBinContent
// + other.fSumWeightsSquared[b] * fUnchangedBinContent * fUnchangedBinContent)
// / (other.fUnchangedBinContent * other.fUnchangedBinContent * other.fUnchangedBinContent * other.fUnchangedBinContent);
// for (size_t b = 0; b < fOverflowSumWeightsSquared.size(); ++b)
// fOverflowSumWeightsSquared[b] = (fOverflowSumWeightsSquared[b] * other.fUnchangedBinContent * other.fUnchangedBinContent
// + other.fOverflowSumWeightsSquared[b] * fUnchangedBinContent * fUnchangedBinContent)
// / (other.fUnchangedBinContent * other.fUnchangedBinContent * other.fUnchangedBinContent * other.fUnchangedBinContent);
}

/// Divide by an other `RHistStatUncertainty` data, assuming same bin configuration.
void DivideBinomial(const RHistStatUncertainty& other) {
Divide(other);
// for (size_t b = 0; b < fSumWeightsSquared.size(); ++b)
// fSumWeightsSquared[b] = std::abs(((1. - 2. * fUnchangedBinContent / other.fUnchangedBinContent) * fSumWeightsSquared[b]
// + fUnchangedBinContent * fUnchangedBinContent * other.fSumWeightsSquared[b] / other.fUnchangedBinContent * other.fUnchangedBinContent)
// / (other.fUnchangedBinContent * other.fUnchangedBinContent));
// for (size_t b = 0; b < fOverflowSumWeightsSquared.size(); ++b)
// fOverflowSumWeightsSquared[b] = std::abs(((1. - 2. * fUnchangedBinContent / other.fUnchangedBinContent) * fOverflowSumWeightsSquared[b]
// + fUnchangedBinContent * fUnchangedBinContent * other.fOverflowSumWeightsSquared[b] / other.fUnchangedBinContent * other.fUnchangedBinContent)
// / (other.fUnchangedBinContent * other.fUnchangedBinContent));
}
};

/** \class RHistDataMomentUncert
Expand Down Expand Up @@ -417,13 +522,41 @@ public:

// FIXME: Add a way to query the inner data

/// Merge with other RHistDataMomentUncert data, assuming same bin configuration.
/// Add an other RHistDataMomentUncert data, assuming same bin configuration.
void Add(const RHistDataMomentUncert& other) {
assert(fMomentXW.size() == other.fMomentXW.size());
assert(fMomentX2W.size() == other.fMomentX2W.size());
for (size_t d = 0; d < DIMENSIONS; ++d) {
fMomentXW[d] += other.fMomentXW[d];
fMomentX2W[d] += other.fMomentX2W[d];
}
}

/// Multiply by a scalar.
template <typename T>
void MultiplyByScalar(const T scalar) {
for (size_t d = 0; d < DIMENSIONS; ++d) {
fMomentXW[d] *= scalar;
fMomentX2W[d] *= scalar;
}
}

/// Divide by an other RHistDataMomentUncert data, assuming same bin configuration.
void Divide(const RHistDataMomentUncert& other) {
// /!\ Should be with the updated weights (after division) of the bin contents
assert(fMomentXW.size() == other.fMomentXW.size());
assert(fMomentX2W.size() == other.fMomentX2W.size());
for (size_t d = 0; d < DIMENSIONS; ++d) {
fMomentXW[d] /= other.fMomentXW[d];
fMomentX2W[d] /= other.fMomentX2W[d];
}
}

/// Divide by an other RHistDataMomentUncert data, assuming same bin configuration.
void DivideBinomial(const RHistDataMomentUncert& other) {
// /!\ Should be with the updated weights (after division) of the bin contents
Divide(other);
}
};

/** \class RHistStatRuntime
Expand Down Expand Up @@ -576,6 +709,41 @@ public:
(void)trigger_base_add{(STAT<DIMENSIONS, PRECISION>::Add(other), 0)...};
}

/// Multiply the current data by a scalar.
template <typename OtherData>
void MultiplyByScalar(const OtherData scalar)
{
// Call `MultiplyByScalar()` on all base classes, using the same tricks as `Fill()`.
using trigger_base_add = int[];
(void)trigger_base_add{(STAT<DIMENSIONS, PRECISION>::MultiplyByScalar(scalar), 0)...};
}

/// Divide current data by some other statistical data.
///
/// The implementation assumes that the other statistics were recorded with
/// the same binning configuration, and that the statistics of `OtherData`
/// are a superset of those recorded by the active `RHistData` instance.
template <typename OtherData>
void Divide(const OtherData &other)
{
// Call `Divide()` on all base classes, using the same tricks as `Fill()`.
using trigger_base_add = int[];
(void)trigger_base_add{(STAT<DIMENSIONS, PRECISION>::Divide(other), 0)...};
}

/// Divide current data by some other statistical data.
///
/// The implementation assumes that the other statistics were recorded with
/// the same binning configuration, and that the statistics of `OtherData`
/// are a superset of those recorded by the active `RHistData` instance.
template <typename OtherData>
void DivideBinomial(const OtherData &other)
{
// Call `DivideBinomial()` on all base classes, using the same tricks as `Fill()`.
using trigger_base_add = int[];
(void)trigger_base_add{(STAT<DIMENSIONS, PRECISION>::DivideBinomial(other), 0)...};
}

/// Whether this provides storage for uncertainties, or whether uncertainties
/// are determined as poisson uncertainty of the content.
static constexpr bool HasBinUncertainty()
Expand Down
12 changes: 0 additions & 12 deletions hist/histv7/inc/ROOT/RHistImpl.hxx
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -42,18 +42,6 @@ using AxisIter_t = std::array<RAxisBase::const_iterator, NDIMS>;
template <int NDIMS>
using AxisIterRange_t = std::array<AxisIter_t<NDIMS>, 2>;

/// Kinds of under- and overflow handling.
enum class EOverflow {
kNoOverflow = 0x0, ///< Exclude under- and overflows
kUnderflow = 0x1, ///< Include underflows
kOverflow = 0x2, ///< Include overflows
kUnderOver = 0x3, ///< Include both under- and overflows
};

inline bool operator&(EOverflow a, EOverflow b)
{
return static_cast<int>(a) & static_cast<int>(b);
}
} // namespace Hist

namespace Detail {
Expand Down
Loading