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

[hageboeck]

The bin computation in TAxis can suffer from floating-point uncertainties, returning bins such that this assertion breaks:

int bin = axis.FindBin(x);
assert(axis.GetBinLowEdge(bin) <= x && x < axis.GetBinUpEdge(bin));

This is of course surprising for users, and sometimes brings data into the overflow, although they fall into the valid range of the axis. Here, two terms are added, which fix these floating-point errors.

On a recent AMD CPU, there is no runtime difference, since the latency of the virtual call overshadows the additional computations for the correction. For 10M calls, running 100 times with perf, we arrive at:

	Uncorrected	Corrected
Instructions	250M	360M
cycles	216M	218M
Run time	0.0495 +- 0.0003 s	0.0493 +- 0.0002 s

Note:
Several attempts at reordering the equations for better precision failed. There were always a few test cases in the included tests that didn't pass, so applying the correction looks to be the best solution.

This PR fixes the issues brought up here:

• https://root-forum.cern.ch/t/floating-point-rounding-error-when-filling-the-histogram/35189
• https://root-forum.cern.ch/t/bug-or-feature-in-ttree-draw/62862/16
• 1703c54
• Bug in TAxis::FindBin (Double_t x) ? #14091

[ferdymercury]

Thanks a lot Stephan!
Some questions from my side:

• Does this supersede [hist] Improve calculation of TAxis::FindBin(x) when x is at the bin edges #14105 ? Maybe it can be closed.
• Does this fix Bug in TAxis::FindBin (Double_t x) ? #14091 ? Maybe it can be linked.
• Does this new implementation affect runtime performance when filling a histogram? If yes, maybe it's better to have a separate function FindBinHighRes or so. One could even think of using long-double.
• Is it sure that the resulting bin is only off by one? I can imagine of (probably very weird cases) where the floating point limited precision might be two?

[hageboeck]

Thanks a lot Stephan! Some questions from my side:

• Does this supersede [hist] Improve calculation of TAxis::FindBin(x) when x is at the bin edges #14105 ? Maybe it can be closed.
• Does this fix Bug in TAxis::FindBin (Double_t x) ? #14091 ? Maybe it can be linked.
• Does this new implementation affect runtime performance when filling a histogram? If yes, maybe it's better to have a separate function FindBinHighRes or so. One could even think of using long-double.
• Is it sure that the resulting bin is only off by one? I can imagine of (probably very weird cases) where the floating point limited precision might be two?

Good questions, I'm checking! 🙂

[github-actions]

Test Results

    19 files      19 suites   4d 3h 1m 0s ⏱️
 3 016 tests  3 013 ✅ 0 💤  3 ❌
55 735 runs  55 708 ✅ 0 💤 27 ❌

For more details on these failures, see this check.

Results for commit 901d041.

♻️ This comment has been updated with latest results.

[hageboeck]

• Does this supersede [hist] Improve calculation of TAxis::FindBin(x) when x is at the bin edges #14105 ? Maybe it can be closed.

It can be closed, indeed.

• Does this fix Bug in TAxis::FindBin (Double_t x) ? #14091 ? Maybe it can be linked.

It fixes the problem described there. I linked the issue.

• Does this new implementation affect runtime performance when filling a histogram? If yes, maybe it's better to have a separate function FindBinHighRes or so. One could even think of using long-double.

It does not affect run time on a modern CPU, even though more instructions are issued. When you run the bare code to find the bin (not as part of a virtual function), you have the following behaviour:

	Uncorrected	Corrected
Instructions	8	19
Cycles	33	31

It doesn't matter how many instructions are issued. What's interesting is how many cycles they need to complete. Although that depends on the CPU in use, you see that the two versions are very close.

The killer argument, however, is that all of this happens inside a virtual function, and that of course has its own latencies. When plugging the codes into a virtual function, and benchmarking 10M calls, we arrive at:

	Uncorrected	Corrected
Instructions	250M	360M
cycles	216M	218M
Run time	0.0495 +- 0.0003 s	0.0493 +- 0.0002 s

So you see that more instructions are issued, but due to the CPU being a superscalar CPU, these additional instructions have no noteworthy effect on the run time, because the correction terms can be evaluated in parallel to the original computation.

This means, that in TAxis, where it is part of a virtual function, we should not see a run time difference.

• Is it sure that the resulting bin is only off by one? I can imagine of (probably very weird cases) where the floating point limited precision might be two?

I did not manage to provoke a more-than-one difference. I didn't set out to prove it mathematically, but I had lots of difficulties crafting a case that made an "off by -1" error.

[ferdymercury]

Awesome analysis, thanks a lot!!
Does this mean that I can close #17689 if you implement the same strategy in that part of the code ?
And not sure if you want to review this alleviation strategy: #17896

[ferdymercury]

Btw, this PR reminds me of this issue: https://its.cern.ch/jira/browse/ROOT-10179

[hageboeck]

Awesome analysis, thanks a lot!! Does this mean that I can close #17689 if you implement the same strategy in that part of the code ?

@ferdymercury I think for this PR, the x == max --> overflow problem remains. So you could probably just use:

- max = currentMax;
+ max = std::nextafter(currentMax, inf);

but that's to be tested.

Btw, this PR reminds me of this issue: https://its.cern.ch/jira/browse/ROOT-10179

Was this the correct link?

[ferdymercury]

Was this the correct link?

It was about the question/fact of whether performance was impacted. Maybe more the subissue https://its.cern.ch/jira/browse/ROOT-10185 which would have been easier for you if rootbench could be triggered from a PR.

 max = std::nextafter(currentMax, inf);
    

I had tried that and it didn't work. I really had to add
xmax = std::max(xmax + 1e-15*(xmax - xmin), std::nextafter(xmax,INFINITY));

For most cases, the 1e-15*max_width is what really counts.
Or do you mean that that should be fine with just nextafter, 'after' your PR gets merged?

Related: #17896