Skip to content

[hist] Improve precision of TAxis::FindFixBin / FindBin. #19033

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

Merged
merged 3 commits into from
Jun 23, 2025
Merged

[hist] Improve precision of TAxis::FindFixBin / FindBin. #19033

Show file tree
Hide file tree
Changes from all commits
Commits
Show all changes
3 commits
Select commit Hold shift + click to select a range
700544e
[NFC] Clean up an outdated comment.
hageboeck Jun 11, 2025
6cb2a8d
[hist] Improve precision of TAxis::FindFixBin / FindBin.
hageboeck Jun 11, 2025
dbe4097
[hist] Add a test for finding bins in TAxis.
hageboeck Jun 11, 2025
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
12 changes: 6 additions & 6 deletions hist/hist/src/TAxis.cxx
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -289,9 +289,6 @@ void TAxis::ExecuteEvent(Int_t event, Int_t px, Int_t py)
/// whereas the overflow bin (`nbins+1`) is for any `x` greater or equal than `fXmax`,
/// as well as for `NaN`. The first regular bin (`1`) is for any `x`
/// greater or equal than `fXmin` and strictly smaller than `fXmin + binwidth`, and so on.
/// @note The bin assignment equation uses doubles, thus rounding errors are
/// expected to appear at the edges. For example: `TAxis(1, -1., 0.).FindBin(-1e-17)`
/// returns the overflow bin (`2`) rather than the theoretically correct bin (`1`).

Int_t TAxis::FindBin(Double_t x)
{
Expand All @@ -314,7 +311,9 @@ Int_t TAxis::FindBin(Double_t x)
return FindFixBin(x);
} else {
if (!fXbins.fN) { //*-* fix bins
bin = 1 + int (fNbins*(x-fXmin)/(fXmax-fXmin) );
const double width = (fXmax-fXmin)/fNbins;
const int approxBin = (x-fXmin) / width; // Might be one unit too small/large due to limited precision of subtraction
bin = 1 + approxBin - (x < fXmin + width*approxBin) + (fXmin + width*(approxBin+1) <= x);
} else { //*-* variable bin sizes
//for (bin =1; x >= fXbins.fArray[bin]; bin++);
bin = 1 + TMath::BinarySearch(fXbins.fN,fXbins.fArray,x);
Expand Down Expand Up @@ -429,9 +428,10 @@ Int_t TAxis::FindFixBin(Double_t x) const
bin = fNbins+1;
} else {
if (!fXbins.fN) { //*-* fix bins
bin = 1 + int (fNbins*(x-fXmin)/(fXmax-fXmin) );
const double width = (fXmax-fXmin)/fNbins;
const int approxBin = (x-fXmin) / width; // Might be one unit too small/large due to limited precision of subtraction
bin = 1 + approxBin - (x < fXmin + width*approxBin) + (fXmin + width*(approxBin+1) <= x);
} else { //*-* variable bin sizes
// for (bin =1; x >= fXbins.fArray[bin]; bin++);
bin = 1 + TMath::BinarySearch(fXbins.fN,fXbins.fArray,x);
}
}
Expand Down
132 changes: 131 additions & 1 deletion hist/hist/test/test_TH1.cxx
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -6,7 +6,7 @@
#include "TH1F.h"
#include "THLimitsFinder.h"


#include <random>
#include <vector>

// StatOverflows TH1
Expand Down Expand Up @@ -149,3 +149,133 @@ TEST(TH1, Normalize)
EXPECT_FLOAT_EQ(v2.Integral("width"), 1.);
EXPECT_FLOAT_EQ(v2.GetMaximum(), 7.9999990);
}

TEST(TAxis, BinComputation_FPAccuracy)
{
// Example from 1703c54
EXPECT_EQ(TAxis(1, -1., 0.).FindBin(-1e-17), 1);

{
// https://root-forum.cern.ch/t/floating-point-rounding-error-when-filling-the-histogram/35189
TAxis axis(128, -0.352, 0.352);
constexpr double x = -0.0220;
EXPECT_EQ(axis.FindBin(x), 61);
EXPECT_EQ(axis.FindFixBin(x), 61);
EXPECT_LE(axis.GetBinLowEdge(61), x);
EXPECT_GT(axis.GetBinUpEdge(61), x);
}

{
// https://github.com/root-project/root/issues/14091
constexpr int nBins = 30;
constexpr double xMin = 3.0, xMax = 6.0;
TAxis ax(nBins, xMin, xMax);

for (Int_t i = 1; i <= ax.GetNbins(); i++) {
EXPECT_EQ(i, ax.FindBin(ax.GetBinLowEdge(i)));
}
}

{
TAxis axis(2000, -1000., 1000.);
EXPECT_EQ(axis.FindFixBin(std::nextafter(-1000., -2000.)), 0);
EXPECT_EQ(axis.FindFixBin(-1000.), 1);
EXPECT_EQ(axis.FindFixBin(std::nextafter(-1000., 0)), 1);

EXPECT_EQ(axis.FindFixBin(-500.00000000001), 500);
EXPECT_EQ(axis.FindFixBin(-500.), 501);
EXPECT_EQ(axis.FindFixBin(-499.9), 501);

EXPECT_EQ(axis.FindFixBin(-1.E-20), 1000);
EXPECT_EQ(axis.FindFixBin(std::nextafter(-0., -1)), 1000);
EXPECT_EQ(axis.FindFixBin(-0.), 1001);
EXPECT_EQ(axis.FindFixBin(0.), 1001);

EXPECT_EQ(axis.FindFixBin(499.9), 1500);
EXPECT_EQ(axis.FindFixBin(500.), 1501);

EXPECT_EQ(axis.FindFixBin(1000. - 1.E-13), 2000);
EXPECT_EQ(axis.FindFixBin(std::nextafter(1000., 0.)), 2000);
EXPECT_EQ(axis.FindFixBin(1000.), 2001);

EXPECT_EQ(axis.FindBin(std::nextafter(-1000., -2000.)), 0);
EXPECT_EQ(axis.FindBin(-1000.), 1);
EXPECT_EQ(axis.FindBin(std::nextafter(-1000., 0)), 1);

EXPECT_EQ(axis.FindBin(-500.00000000001), 500);
EXPECT_EQ(axis.FindBin(-500.), 501);
EXPECT_EQ(axis.FindBin(-499.9), 501);

EXPECT_EQ(axis.FindBin(-1.E-20), 1000);
EXPECT_EQ(axis.FindBin(std::nextafter(-0., -1)), 1000);
EXPECT_EQ(axis.FindBin(-0.), 1001);
EXPECT_EQ(axis.FindBin(0.), 1001);

EXPECT_EQ(axis.FindBin(499.9), 1500);
EXPECT_EQ(axis.FindBin(500.), 1501);

EXPECT_EQ(axis.FindBin(1000. - 1.E-13), 2000);
EXPECT_EQ(axis.FindBin(std::nextafter(1000., 0.)), 2000);
EXPECT_EQ(axis.FindBin(1000.), 2001);
}

for (const auto &[low, high] : std::initializer_list<std::pair<long double, long double>>{{-10654.1l, 32165.l},
{-1.0656E23l, -20654.l},
{1.1234E4l, 4.5678E20l},
{1.E-60l, 1.E-20l},
{-1.E-20l, -1.E-60l}}) {
constexpr int N = 100;
const double width = (high - low) / N;
std::mt19937 gen;
std::uniform_real_distribution dist{low - width, high + width};
TAxis axis(N, low, high);
for (unsigned int i = 0; i < 100000; ++i) {
long double x = dist(gen);
if (i == 0) {
x = std::nextafter(double(high), double(low));
}
const auto bin = axis.FindFixBin(x);
EXPECT_EQ(bin, axis.FindBin(x));

if (x < double(low)) {
EXPECT_EQ(bin, 0);
} else if (x >= double(high)) {
EXPECT_EQ(bin, N + 1);
}

EXPECT_LE(axis.GetBinLowEdge(bin), x);
EXPECT_LT(x, axis.GetBinUpEdge(bin));
}
}
}

// The tests below are taken from https://github.com/root-project/root/pull/14105
TEST(TAxis, FindBinExact)
{
// Test the case where bin edges are exactly represented as floating points
TAxis ax(88, 1010, 1098);
for (int i = 1; i <= ax.GetNbins(); i++) {
double x = ax.GetBinLowEdge(i);
EXPECT_EQ(i, ax.FindBin(x));
EXPECT_EQ(i, ax.FindFixBin(x));
x = ax.GetBinUpEdge(i);
EXPECT_EQ(i + 1, ax.FindBin(x));
EXPECT_EQ(i + 1, ax.FindFixBin(x));
x -= x * std::numeric_limits<double>::epsilon();
EXPECT_EQ(i, ax.FindBin(x));
}
}
TEST(TAxis, FindBinApprox)
{
TAxis ax(90, 0., 10.);
for (int i = 1; i <= ax.GetNbins(); i++) {
double x = ax.GetBinLowEdge(i);
EXPECT_EQ(i, ax.FindBin(x));
EXPECT_EQ(i, ax.FindFixBin(x));
x = ax.GetBinUpEdge(i);
EXPECT_EQ(i + 1, ax.FindBin(x));
EXPECT_EQ(i + 1, ax.FindFixBin(x));
x -= x * std::numeric_limits<double>::epsilon();
EXPECT_EQ(i, ax.FindBin(x));
}
}
Loading