Ticket #8085 (closed: fixed)
Problem with Integration Boundaries
Reported by: | Owen Arnold | Owned by: | Federico M Pouzols |
---|---|---|---|
Priority: | critical | Milestone: | Release 3.3 |
Component: | Framework | Keywords: | NewStarter |
Cc: | Blocked By: | ||
Blocking: | Tester: | Dan Nixon |
Description
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.
Change History
comment:2 Changed 6 years ago by Nick Draper
- Keywords NewStarter added
- Owner changed from Anyone to Nick Draper
comment:3 Changed 6 years ago by Nick Draper
- Owner changed from Nick Draper to Federico M Pouzols
- Milestone changed from Backlog to Release 3.3
comment:4 Changed 6 years ago by Federico Montesino Pouzols
- Status changed from assigned to inprogress
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
comment:5 Changed 6 years ago by 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
comment:6 Changed 6 years ago by Federico Montesino Pouzols
find RangeLower with std::lower_bound, neat, re #8085
This also fixes some issues in ChopDataTest
Changeset: 8d091cb521f0a8ccbdce743327cf9083611aaca3
comment:7 Changed 6 years ago by Federico Montesino Pouzols
New tests for Integration alg, (ab)using real numbers, re #8085
Changeset: 7b132fa8b60562b3578e7ef1870b56b454617736
comment:8 Changed 6 years ago by Federico Montesino Pouzols
constness fix, re #8085
Changeset: d3191486ec677c8dcfed59b2396fedaf267a796e
comment:9 Changed 6 years ago by Federico M Pouzols
- Status changed from inprogress to verify
- Resolution set to fixed
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).
comment:10 Changed 6 years ago by Dan Nixon
- Status changed from verify to verifying
- Tester set to Dan Nixon
comment:11 Changed 6 years ago by Dan Nixon
- Status changed from verifying to closed
Merge remote-tracking branch 'origin/bugfix/8085_integration_boundaries'
Full changeset: 711f365241cdb0a1c0e24f84628ec95ba21d8607
comment:12 Changed 5 years ago by Stuart Campbell
This ticket has been transferred to github issue 8930
Bulk move to assigned at the introduction of the triage step