Problem with Integration Boundaries

[Owen Arnold]

try this:

CreateWorkspace(OutputWorkspace='ws2',DataX='-1,-0.8,-0.6,-0.4,-0.2,-2.22045e-16,0.2,0.4,0.6,0.8,1',DataY='0,0,0,2,2,2,2,0,0,0',DataE='0,0,0,0,0,0,0,0,0,0',UnitX='1/q')
Rebin(InputWorkspace='ws2',OutputWorkspace='rhs_rebinned',Params='0.2')
Integration(InputWorkspace='rhs_rebinned',OutputWorkspace='rhsOverlapIntegrated',RangeLower=-0.4,RangeUpper=-0.2)

The output workspace (rhsOverlapIntegrated) should have a single value of 2, because I specified boundaries that would exactly result in a single bin being captured and an integrated value of 2.0 not 0.0. But it doesn't do this, so It's not checking the boundaries properly!

As possible further evidence for this, look at the history generated from rhsOverlapIntegrated and you'll see that the RangeLower and RangeUpper have values being reported that are not exactly what I specified. Could this be the cause of the issue.

[Nick Draper]

Bulk move to assigned at the introduction of the triage step

[Federico Montesino Pouzols]

Fix issue, based on find_if, , re #8085

This is now more robust to numerical errors in bin limits, it passes all tests, including new ones.

Changeset: a12312cc6b31b2296817e7d377d64a6c3b23ef04

[Federico Montesino Pouzols]

RangeUpper: use std::upper_bound rather than find_if, re #8085

upper_bound will be faster if any, and find_if is deprecated

Changeset: 91dff39e3e6c4458d828111739ddcfb3bca89b68

[Federico Montesino Pouzols]

find RangeLower with std::lower_bound, neat, re #8085

This also fixes some issues in ChopDataTest

Changeset: 8d091cb521f0a8ccbdce743327cf9083611aaca3

[Federico Montesino Pouzols]

New tests for Integration alg, (ab)using real numbers, re #8085

Changeset: 7b132fa8b60562b3578e7ef1870b56b454617736

[Federico Montesino Pouzols]

constness fix, re #8085

Changeset: d3191486ec677c8dcfed59b2396fedaf267a796e

[Federico M Pouzols]

This is now working and passing all tests. I've extended the IntegrationTest unit test with new tests based on the example given in the ticket description. The tests were already quite exhaustive but what was keeping this bug/issue hidden is that most tests define the bin boundaries as integers, and integers can be represented exactly in double values.

Things get worse when the boundaries are real numbers with some precision and epsilon/error. I've changed the lines that search for LowerRange and UpperRange, the algorithms used now should be more robust, and they also use a "tolerant" comparator (to an error of 1 double epsilon).

To test: for a start, you can test the example python commands given in the ticket, now it should produce an output of 2 (you can check this in the Y values of the Integration output workspace, or 'rhs_rebinned' in the example). You can also test variations of it, such as integrating between -0.2 and 0.2 (should give 4), etc. There are also new tests included in IntegrationTest. The most relevant unit tests here are AlgorithmsTest_IntegrationTest and AlgorithmsTest_ChopDataTest:

ctest -V -R AlgorithmsTest_IntegrationTest
ctest -V -R AlgorithmsTest_ChopDataTest

---

---

And here is a bit more verbose explanation of what was going on. This is a delicate issue. I think this particular ticket may be closed but the more general issue it raises maybe not. We may need some discussion about the more general issue: what do we do with floating point precision issues in algorithms that iterate through / depend on boundaries or thresholds.

In this particular case, with Integration, if you do

 Integration(InputWorkspace='rhs_rebinned',OutputWorkspace='test_rhsOverlapIntegrated',RangeLower=-0.4,RangeUpper=-0.2)

In the algorithm history of the test_rhsOverlapIntegrated workspace you'll see in the properties of the integration algorithm that:

 RangeLower: -0.40000000000000002
 RangeUpper: -0.20000000000000001

And the output is 0 (wrong)

Now, if you try:

 Integration(InputWorkspace='rhs_rebinned',OutputWorkspace='test_rhsOverlapIntegrated',RangeLower=-0.4000000001,RangeUpper=-0.2)

(or similar, just put a few 0s between the 4 and the 1, but less than 15 0s)

The output is 2 (correct).

Changing the upper boundary, from -0.2 to -0.19999999, etc. has no effect.

So at first it seems there is an issue with the precision of C double real numbers, as hinted in the ticket description.

When you debug the algorithm, what you find is that RangeLower and RangeUpper are fine (within the double type's epsilon). But lowit (see Code/Mantid/Framework/Algorithms/src/Integration.cpp: Integration:exec()) was being searched in a way that is sensitive to this type of numerical issues (as opposed to highit).

highit was found as the first bin boundary that is greater than 'RangeUpper', then decreasing the iterator to the previous bin boundary. However lowit uses std::lower_bound which finds the first element that is not less than 'RangeLower'. RangeUpper and RangeLower are now found using std::upper_bound and std::upper_lower with a custom comparator that sees values differing by just one double epsilon as equal (or not less than).

[Dan Nixon]

Merge remote-tracking branch 'origin/bugfix/8085_integration_boundaries'

Full changeset: ​711f365241cdb0a1c0e24f84628ec95ba21d8607

[Stuart Campbell]

This ticket has been transferred to github issue ​8930