diff --git a/Framework/Algorithms/inc/MantidAlgorithms/CompareWorkspaces.h b/Framework/Algorithms/inc/MantidAlgorithms/CompareWorkspaces.h index d41c1960b9bd..c17c71713836 100644 --- a/Framework/Algorithms/inc/MantidAlgorithms/CompareWorkspaces.h +++ b/Framework/Algorithms/inc/MantidAlgorithms/CompareWorkspaces.h @@ -76,8 +76,8 @@ class MANTID_ALGORITHMS_DLL CompareWorkspaces final : public API::Algorithm { "testing process."; } - static bool withinAbsoluteTolerance(double x1, double x2, double atol); - static bool withinRelativeTolerance(double x1, double x2, double rtol); + static bool withinAbsoluteTolerance(double x1, double x2, double atol, bool const nanEqual = false); + static bool withinRelativeTolerance(double x1, double x2, double rtol, bool const nanEqual = false); private: /// Initialise algorithm diff --git a/Framework/Algorithms/src/CompareWorkspaces.cpp b/Framework/Algorithms/src/CompareWorkspaces.cpp index 9119d8064bdc..56324ed3c343 100644 --- a/Framework/Algorithms/src/CompareWorkspaces.cpp +++ b/Framework/Algorithms/src/CompareWorkspaces.cpp @@ -22,6 +22,7 @@ #include "MantidGeometry/Crystal/IPeak.h" #include "MantidGeometry/Instrument/ComponentInfo.h" #include "MantidGeometry/Instrument/DetectorInfo.h" +#include "MantidKernel/FloatingPointComparison.h" #include "MantidKernel/Unit.h" namespace Mantid::Algorithms { @@ -71,23 +72,18 @@ int compareEventLists(Kernel::Logger &logger, const EventList &el1, const EventL const auto &e1 = events1[i]; const auto &e2 = events2[i]; - bool diffpulse = false; - bool difftof = false; - bool diffweight = false; - if (std::abs(e1.pulseTime().totalNanoseconds() - e2.pulseTime().totalNanoseconds()) > tolPulse) { - diffpulse = true; - ++numdiffpulse; - } - if (std::abs(e1.tof() - e2.tof()) > tolTof) { - difftof = true; - ++numdifftof; - } + bool diffpulse = + !withinAbsoluteDifference(e1.pulseTime().totalNanoseconds(), e2.pulseTime().totalNanoseconds(), tolPulse); + bool difftof = !withinAbsoluteDifference(e1.tof(), e2.tof(), tolTof); + bool diffweight = !withinAbsoluteDifference(e1.weight(), e2.weight(), tolWeight); if (diffpulse && difftof) - ++numdiffboth; - if (std::abs(e1.weight() - e2.weight()) > tolWeight) { - diffweight = true; - ++numdiffweight; - } + numdiffboth++; + if (diffpulse) + numdiffpulse++; + if (difftof) + numdifftof++; + if (diffweight) + numdiffweight++; bool same = (!diffpulse) && (!difftof) && (!diffweight); if (!same) { @@ -148,6 +144,8 @@ void CompareWorkspaces::init() { "Very often such logs are huge so making it true should be " "the last option."); + declareProperty("NaNsEqual", false, "Whether NaN values should compare as equal with other NaN values."); + declareProperty("NumberMismatchedSpectraToPrint", 1, "Number of mismatched spectra from lowest to be listed. "); declareProperty("DetailedPrintIndex", EMPTY_INT(), "Mismatched spectra that will be printed out in details. "); @@ -172,13 +170,14 @@ void CompareWorkspaces::exec() { m_parallelComparison = false; double const tolerance = getProperty("Tolerance"); + bool const nanEqual = getProperty("NaNsEqual"); if (getProperty("ToleranceRelErr")) { - this->m_compare = [tolerance](double const x1, double const x2) -> bool { - return CompareWorkspaces::withinRelativeTolerance(x1, x2, tolerance); + this->m_compare = [tolerance, nanEqual](double const x1, double const x2) -> bool { + return CompareWorkspaces::withinRelativeTolerance(x1, x2, tolerance, nanEqual); }; } else { - this->m_compare = [tolerance](double const x1, double const x2) -> bool { - return CompareWorkspaces::withinAbsoluteTolerance(x1, x2, tolerance); + this->m_compare = [tolerance, nanEqual](double const x1, double const x2) -> bool { + return CompareWorkspaces::withinAbsoluteTolerance(x1, x2, tolerance, nanEqual); }; } @@ -1049,6 +1048,7 @@ void CompareWorkspaces::doPeaksComparison(PeaksWorkspace_sptr tws1, PeaksWorkspa } const bool isRelErr = getProperty("ToleranceRelErr"); + const bool checkAllData = getProperty("CheckAllData"); for (int i = 0; i < tws1->getNumberPeaks(); i++) { const Peak &peak1 = tws1->getPeak(i); const Peak &peak2 = tws2->getPeak(i); @@ -1127,7 +1127,8 @@ void CompareWorkspaces::doPeaksComparison(PeaksWorkspace_sptr tws1, PeaksWorkspa << "value1 = " << s1 << "\n" << "value2 = " << s2 << "\n"; recordMismatch("Data mismatch"); - return; + if (!checkAllData) + return; } } } @@ -1163,6 +1164,8 @@ void CompareWorkspaces::doLeanElasticPeaksComparison(const LeanElasticPeaksWorks const double tolerance = getProperty("Tolerance"); const bool isRelErr = getProperty("ToleranceRelErr"); + const bool checkAllData = getProperty("CheckAllData"); + const bool nanEqual = getProperty("NaNsEqual"); for (int peakIndex = 0; peakIndex < ipws1->getNumberPeaks(); peakIndex++) { for (size_t j = 0; j < ipws1->columnCount(); j++) { std::shared_ptr col = ipws1->getColumn(j); @@ -1229,10 +1232,10 @@ void CompareWorkspaces::doLeanElasticPeaksComparison(const LeanElasticPeaksWorks // bool mismatch = !m_compare(s1, s2) // can replace this if/else, and isRelErr and tolerance can be deleted if (isRelErr && name != "QLab" && name != "QSample") { - if (!withinRelativeTolerance(s1, s2, tolerance)) { + if (!withinRelativeTolerance(s1, s2, tolerance, nanEqual)) { mismatch = true; } - } else if (!withinAbsoluteTolerance(s1, s2, tolerance)) { + } else if (!withinAbsoluteTolerance(s1, s2, tolerance, nanEqual)) { mismatch = true; } if (mismatch) { @@ -1242,7 +1245,8 @@ void CompareWorkspaces::doLeanElasticPeaksComparison(const LeanElasticPeaksWorks << "value1 = " << s1 << "\n" << "value2 = " << s2 << "\n"; recordMismatch("Data mismatch"); - return; + if (!checkAllData) + return; } } } @@ -1356,12 +1360,15 @@ this error is within the limits requested. @param x1 -- first value to check difference @param x2 -- second value to check difference @param atol -- the tolerance of the comparison. Must be nonnegative +@param nanEqual -- whether two NaNs compare as equal @returns true if absolute difference is within the tolerance; false otherwise */ -bool CompareWorkspaces::withinAbsoluteTolerance(double const x1, double const x2, double const atol) { - // NOTE !(|x1-x2| > atol) is not the same as |x1-x2| <= atol - return !(std::abs(x1 - x2) > atol); +bool CompareWorkspaces::withinAbsoluteTolerance(double const x1, double const x2, double const atol, + bool const nanEqual) { + if (nanEqual && std::isnan(x1) && std::isnan(x2)) + return true; + return Kernel::withinAbsoluteDifference(x1, x2, atol); } //------------------------------------------------------------------------------------------------ @@ -1371,24 +1378,15 @@ this error is within the limits requested. @param x1 -- first value to check difference @param x2 -- second value to check difference @param rtol -- the tolerance of the comparison. Must be nonnegative +@param nanEqual -- whether two NaNs compare as equal @returns true if relative difference is within the tolerance; false otherwise @returns true if error or false if the relative value is within the limits requested */ -bool CompareWorkspaces::withinRelativeTolerance(double const x1, double const x2, double const rtol) { - // calculate difference - double const num = std::abs(x1 - x2); - // return early if the values are equal - if (num == 0.0) +bool CompareWorkspaces::withinRelativeTolerance(double const x1, double const x2, double const rtol, + bool const nanEqual) { + if (nanEqual && std::isnan(x1) && std::isnan(x2)) return true; - // create the average magnitude for comparison - double const den = 0.5 * (std::abs(x1) + std::abs(x2)); - // return early, possibly avoids a multiplication - // NOTE if den<1, then divsion will only make num larger - // NOTE if den<1 but num<=rtol, we cannot conclude anything - if (den <= 1.0 && num > rtol) - return false; - // NOTE !(num > rtol*den) is not the same as (num <= rtol*den) - return !(num > (rtol * den)); + return Kernel::withinRelativeDifference(x1, x2, rtol); } } // namespace Mantid::Algorithms diff --git a/Framework/Algorithms/src/DetectorEfficiencyCor.cpp b/Framework/Algorithms/src/DetectorEfficiencyCor.cpp index c21580d145d2..ae7187fe744d 100644 --- a/Framework/Algorithms/src/DetectorEfficiencyCor.cpp +++ b/Framework/Algorithms/src/DetectorEfficiencyCor.cpp @@ -123,7 +123,7 @@ void DetectorEfficiencyCor::exec() { int64_t numHists = m_inputWS->getNumberHistograms(); auto numHists_d = static_cast(numHists); const auto progStep = static_cast(ceil(numHists_d / 100.0)); - auto &spectrumInfo = m_inputWS->spectrumInfo(); + auto const &spectrumInfo = m_inputWS->spectrumInfo(); PARALLEL_FOR_IF(Kernel::threadSafe(*m_inputWS, *m_outputWS)) for (int64_t i = 0; i < numHists; ++i) { diff --git a/Framework/Algorithms/test/CompareWorkspacesTest.h b/Framework/Algorithms/test/CompareWorkspacesTest.h index d97b3b71fa8b..44fc3da12869 100644 --- a/Framework/Algorithms/test/CompareWorkspacesTest.h +++ b/Framework/Algorithms/test/CompareWorkspacesTest.h @@ -1193,6 +1193,60 @@ class CompareWorkspacesTest : public CxxTest::TestSuite { TS_ASSERT_EQUALS(alg.getPropertyValue("Result"), PROPERTY_VALUE_TRUE); } + void test_equal_tableworkspaces_match() { + std::string const col_type("double"), col_name("aColumn"); + std::vector col_values{1.0, 2.0, 3.0}; + // create the table workspaces + Mantid::API::ITableWorkspace_sptr table1 = WorkspaceFactory::Instance().createTable(); + table1->addColumn(col_type, col_name); + for (double val : col_values) { + TableRow newrow = table1->appendRow(); + newrow << val; + } + Mantid::API::ITableWorkspace_sptr table2 = WorkspaceFactory::Instance().createTable(); + table2->addColumn(col_type, col_name); + for (double val : col_values) { + TableRow newrow = table2->appendRow(); + newrow << val; + } + + Mantid::Algorithms::CompareWorkspaces alg; + alg.initialize(); + TS_ASSERT_THROWS_NOTHING(alg.setProperty("Workspace1", table1)); + TS_ASSERT_THROWS_NOTHING(alg.setProperty("Workspace2", table2)); + TS_ASSERT(alg.execute()); + TS_ASSERT_EQUALS(alg.getPropertyValue("Result"), PROPERTY_VALUE_TRUE); + } + + void test_tableworkspace_NaNs_fails() { + std::string const col_type("double"), col_name("aColumn"); + std::vector col_values1{1.0, 2.0, 3.0}; + std::vector col_values2{1.0, 2.0, std::numeric_limits::quiet_NaN()}; + // create the table workspaces + Mantid::API::ITableWorkspace_sptr table1 = WorkspaceFactory::Instance().createTable(); + table1->addColumn(col_type, col_name); + for (double val : col_values1) { + TableRow newrow = table1->appendRow(); + newrow << val; + } + Mantid::API::ITableWorkspace_sptr table2 = WorkspaceFactory::Instance().createTable(); + table2->addColumn(col_type, col_name); + for (double val : col_values2) { + TableRow newrow = table2->appendRow(); + newrow << val; + } + + Mantid::Algorithms::CompareWorkspaces alg; + alg.initialize(); + TS_ASSERT_THROWS_NOTHING(alg.setProperty("Workspace1", table1)); + TS_ASSERT_THROWS_NOTHING(alg.setProperty("Workspace2", table2)); + TS_ASSERT(alg.execute()); + TS_ASSERT_EQUALS(alg.getPropertyValue("Result"), PROPERTY_VALUE_FALSE); + + ITableWorkspace_sptr table = AnalysisDataService::Instance().retrieveWS("compare_msgs"); + TS_ASSERT_EQUALS(table->cell(0, 0), "Table data mismatch"); + } + void test_tableworkspace_different_column_names_fails() { auto table1 = setupTableWorkspace(); table1->getColumn(5)->setName("SomethingElse"); diff --git a/Framework/Algorithms/test/Q1DWeightedTest.h b/Framework/Algorithms/test/Q1DWeightedTest.h index 1d845f05c640..7dcfa5b0b008 100644 --- a/Framework/Algorithms/test/Q1DWeightedTest.h +++ b/Framework/Algorithms/test/Q1DWeightedTest.h @@ -29,6 +29,11 @@ using Mantid::DataHandling::LoadNexusProcessed; using Mantid::DataHandling::LoadSpice2D; using Mantid::DataHandling::MoveInstrumentComponent; +namespace { +bool const USE_NANS_EQUAL(true); +bool const USE_NANS_NOT_EQUAL(false); +} // namespace + class Q1DWeightedTest : public CxxTest::TestSuite { public: void testName() { TS_ASSERT_EQUALS(radial_average.name(), "Q1DWeighted") } @@ -277,6 +282,11 @@ class Q1DWeightedTest : public CxxTest::TestSuite { compareWorkspaces(refWedges, outputWedges); } + /** + * The result and the expected value used in this test are two matrix + * workspaces with NaNs in all y-values and in all e-values. + * Is this even a meaningful test? + */ void testShapeTableResultsAsymm() { // test the results computed by the table shape method against those from // the usual method when asymmetricWedges is set to true @@ -295,7 +305,7 @@ class Q1DWeightedTest : public CxxTest::TestSuite { populateAlgorithm(refWS, refWedges, false, true, 4); TS_ASSERT_THROWS_NOTHING(radial_average.execute()) - compareWorkspaces(refWedges, outputWedges); + compareWorkspaces(refWedges, outputWedges, USE_NANS_EQUAL); } void testShapeCorrectOrder() { @@ -509,7 +519,7 @@ class Q1DWeightedTest : public CxxTest::TestSuite { } } - void compareWorkspaces(std::string refWS, std::string toCompare) { + void compareWorkspaces(std::string refWS, std::string toCompare, bool nansEqual = USE_NANS_NOT_EQUAL) { WorkspaceGroup_sptr result; TS_ASSERT_THROWS_NOTHING( result = std::dynamic_pointer_cast(AnalysisDataService::Instance().retrieve(toCompare))) @@ -529,6 +539,7 @@ class Q1DWeightedTest : public CxxTest::TestSuite { comparison.setPropertyValue("CheckAllData", "1"); comparison.setPropertyValue("CheckType", "1"); comparison.setPropertyValue("ToleranceRelErr", "1"); + comparison.setProperty("NaNsEqual", nansEqual); TS_ASSERT(comparison.execute()) TS_ASSERT(comparison.isExecuted()) TS_ASSERT_EQUALS(comparison.getPropertyValue("Result"), "1"); diff --git a/Framework/DataObjects/inc/MantidDataObjects/TableColumn.h b/Framework/DataObjects/inc/MantidDataObjects/TableColumn.h index 752601e0b5e6..52a3c70210bf 100644 --- a/Framework/DataObjects/inc/MantidDataObjects/TableColumn.h +++ b/Framework/DataObjects/inc/MantidDataObjects/TableColumn.h @@ -14,6 +14,7 @@ #include #include "MantidAPI/Column.h" +#include "MantidKernel/FloatingPointComparison.h" #include "MantidKernel/V3D.h" namespace Mantid { @@ -64,10 +65,10 @@ template class TableColumn : public API::Column { std::string name = std::string(typeid(Type).name()); if ((name.find('i') != std::string::npos) || (name.find('l') != std::string::npos) || (name.find('x') != std::string::npos)) { - if (length == 4) { + if (length == 4) { // cppcheck-suppress knownConditionTrueFalse this->m_type = "int"; } - if (length == 8) { + if (length == 8) { // cppcheck-suppress knownConditionTrueFalse this->m_type = "int64"; } } @@ -78,10 +79,10 @@ template class TableColumn : public API::Column { this->m_type = "double"; } if (name.find('u') != std::string::npos) { - if (length == 4) { + if (length == 4) { // cppcheck-suppress knownConditionTrueFalse this->m_type = "uint32_t"; } - if (length == 8) { + if (length == 8) { // cppcheck-suppress knownConditionTrueFalse this->m_type = "uint64_t"; } } @@ -221,7 +222,7 @@ template class TableColumn : public API::Column { // helper function template for equality bool compareVectors(const std::vector &newVector, double tolerance) const { for (size_t i = 0; i < m_data.size(); i++) { - if (fabs((double)m_data[i] - (double)newVector[i]) > tolerance) { + if (!Kernel::withinAbsoluteDifference(m_data[i], newVector[i], tolerance)) { return false; } } @@ -231,17 +232,14 @@ template class TableColumn : public API::Column { // helper function template for equality with relative error bool compareVectorsRelError(const std::vector &newVector, double tolerance) const { for (size_t i = 0; i < m_data.size(); i++) { - double num = fabs((double)m_data[i] - (double)newVector[i]); - double den = (fabs((double)m_data[i]) + fabs((double)newVector[i])) / 2; - if (den < tolerance && num > tolerance) { - return false; - } else if (num / den > tolerance) { + if (!Kernel::withinRelativeDifference(m_data[i], newVector[i], tolerance)) { return false; } } return true; } }; + /// Template specialisation for long64 template <> inline bool TableColumn::compareVectors(const std::vector &newVector, double tolerance) const { diff --git a/Framework/Kernel/inc/MantidKernel/FloatingPointComparison.h b/Framework/Kernel/inc/MantidKernel/FloatingPointComparison.h index eb4ee09791fc..b4e562dda304 100644 --- a/Framework/Kernel/inc/MantidKernel/FloatingPointComparison.h +++ b/Framework/Kernel/inc/MantidKernel/FloatingPointComparison.h @@ -10,6 +10,7 @@ // Includes //----------------------------------------------------------------------------- #include "MantidKernel/DllConfig.h" +#include "MantidKernel/V3D.h" namespace Mantid { namespace Kernel { @@ -24,8 +25,10 @@ template MANTID_KERNEL_DLL T absoluteDifference(T const x, T const /// Calculate relative difference between x, y template MANTID_KERNEL_DLL T relativeDifference(T const x, T const y); /// Test whether x, y are within absolute tolerance tol -template MANTID_KERNEL_DLL bool withinAbsoluteDifference(T const x, T const y, T const tolerance); +template +MANTID_KERNEL_DLL bool withinAbsoluteDifference(T const x, T const y, S const tolerance); /// Test whether x, y are within relative tolerance tol -template MANTID_KERNEL_DLL bool withinRelativeDifference(T const x, T const y, T const tolerance); +template +MANTID_KERNEL_DLL bool withinRelativeDifference(T const x, T const y, S const tolerance); } // namespace Kernel } // namespace Mantid diff --git a/Framework/Kernel/src/FloatingPointComparison.cpp b/Framework/Kernel/src/FloatingPointComparison.cpp index f2c15a5492e1..bd694413b320 100644 --- a/Framework/Kernel/src/FloatingPointComparison.cpp +++ b/Framework/Kernel/src/FloatingPointComparison.cpp @@ -9,6 +9,7 @@ //----------------------------------------------------------------------------- #include "MantidKernel/FloatingPointComparison.h" #include "MantidKernel/System.h" +#include "MantidKernel/V3D.h" #include #include @@ -20,10 +21,26 @@ namespace Mantid::Kernel { * std::numeric_limits::epsilon precision * @param x :: LHS comparator * @param y :: RHS comparator - * @returns True if the numbers are considered equal within the given tolerance, - * false otherwise. False if any value is NaN. + * @returns True if the numbers are considered equal within machine precision. + * Machine precision is a 1 at the least significant bit scaled to the same power as the numbers compared. + * E.g. for 1, it is usually 1x2^{-52} + * E.g. for 1x2^100, it is usually 1x2^{-52} x 1x2^100 = 1x2^{48} + * False if any value is NaN. False if comparing opposite infinities. */ -template bool equals(T const x, T const y) { return std::abs(x - y) <= std::numeric_limits::epsilon(); } +template bool equals(T const x, T const y) { + // handle infinities + if (std::isinf(x) && std::isinf(y)) { + // if x,y both +inf, x-y=NaN; if x,y both -inf, x-y=NaN; else not an NaN + return std::isnan(x - y); + } else { + // produce a scale for comparison + // in general can use either value, but use the second to work better with near differences, + // which call as equals(difference, tolerance), since tolerance will more often be finite + int const exp = y < std::numeric_limits::min() ? std::numeric_limits::min_exponent - 1 : std::ilogb(y); + // compare to within machine epsilon + return std::abs(x - y) <= std::ldexp(std::numeric_limits::epsilon(), exp); + } +} /** * Compare two floating-point numbers as to whether they satisfy x<=y within @@ -68,10 +85,10 @@ template T relativeDifference(T const x, T const y) { T const num = absoluteDifference(x, y); if (num <= std::numeric_limits::epsilon()) { // if |x-y| == 0.0 (within machine tolerance), relative difference is zero - return 0.0; + return static_cast(0); } else { // otherwise we have to calculate the denominator - T const denom = static_cast(0.5 * (std::abs(x) + std::abs(y))); + T const denom = static_cast((std::abs(x) + std::abs(y)) / static_cast(2)); // NOTE if we made it this far, at least one of x or y is nonzero, so denom will be nonzero return num / denom; } @@ -86,8 +103,12 @@ template T relativeDifference(T const x, T const y) { * @returns True if the numbers are considered equal within the given tolerance, * false otherwise. False if either value is NaN. */ -template bool withinAbsoluteDifference(T const x, T const y, T const tolerance) { - return ltEquals(absoluteDifference(x, y), tolerance); +template bool withinAbsoluteDifference(T const x, T const y, S const tolerance) { + // handle the case of infinities + if (std::isinf(x) && std::isinf(y)) + // if both are +inf, return true; if both -inf, return true; else false + return isnan(static_cast(x - y)); + return ltEquals(static_cast(absoluteDifference(x, y)), tolerance); } /** @@ -99,21 +120,21 @@ template bool withinAbsoluteDifference(T const x, T const y, T cons * @returns True if the numbers are considered equal within the given tolerance, * false otherwise. False if either value is NaN. */ -template bool withinRelativeDifference(T const x, T const y, T const tolerance) { - // handle the case of NaNs - if (std::isnan(x) || std::isnan(y)) { - return false; - } - T const num = absoluteDifference(x, y); - if (!(num > std::numeric_limits::epsilon())) { +template bool withinRelativeDifference(T const x, T const y, S const tolerance) { + // handles the case of infinities + if (std::isinf(x) && std::isinf(y)) + // if both are +inf, return true; if both -inf, return true; else false + return isnan(static_cast(x - y)); + S const num = static_cast(absoluteDifference(x, y)); + if (num <= std::numeric_limits::epsilon()) { // if |x-y| == 0.0 (within machine tolerance), this test passes return true; } else { // otherwise we have to calculate the denominator - T const denom = static_cast(0.5 * (std::abs(x) + std::abs(y))); + S const denom = static_cast(std::abs(x) + std::abs(y)) / static_cast(2); // if denom <= 1, then |x-y| > tol implies |x-y|/denom > tol, can return early // NOTE can only return early if BOTH denom > tol AND |x-y| > tol. - if (denom <= 1. && !ltEquals(num, tolerance)) { + if (denom <= static_cast(1) && !ltEquals(num, tolerance)) { return false; } else { // avoid division for improved performance @@ -130,16 +151,59 @@ template DLLExport bool ltEquals(double const, double const); template DLLExport bool ltEquals(float const, float const); template DLLExport bool gtEquals(double const, double const); template DLLExport bool gtEquals(float const, float const); -// +// difference methods template DLLExport double absoluteDifference(double const, double const); template DLLExport float absoluteDifference(float const, float const); template DLLExport double relativeDifference(double const, double const); template DLLExport float relativeDifference(float const, float const); -// +// within difference methods template DLLExport bool withinAbsoluteDifference(double const, double const, double const); template DLLExport bool withinAbsoluteDifference(float const, float const, float const); +template DLLExport bool withinAbsoluteDifference(long const, long const, long const); +template DLLExport bool withinAbsoluteDifference(long long const, long long const, long long const); template DLLExport bool withinRelativeDifference(double const, double const, double const); template DLLExport bool withinRelativeDifference(float const, float const, float const); +template DLLExport bool withinRelativeDifference(long const, long const, long const); +template DLLExport bool withinRelativeDifference(long long const, long long const, long long const); +// instantiations where the objects are anything and tolerance is a double +template DLLExport bool withinAbsoluteDifference(float const, float const, double const); +template DLLExport bool withinAbsoluteDifference(int const, int const, double const); +template DLLExport bool withinRelativeDifference(float const, float const, double const); +template DLLExport bool withinRelativeDifference(int const, int const, double const); ///@endcond +/// Template specialization for long int +/// these prevent compiler warnings about using isinf and isnan on longs + +template <> DLLExport bool equals(long const x, long const y) { return x == y; } + +template <> DLLExport bool withinAbsoluteDifference(long const x, long const y, long const tolerance) { + return ltEquals(std::labs(x - y), tolerance); +} + +template <> DLLExport bool withinRelativeDifference(long const x, long const y, long const tolerance) { + long const num = std::labs(x - y); + if (num == 0) { + return true; + } else { + long const denom = (std::labs(x) + std::labs(y)) / 2; + return num <= static_cast(denom * tolerance); + } +} + +/// Template specialization for unsigned int + +template <> +DLLExport bool withinAbsoluteDifference(unsigned int const x, unsigned int const y, + double const tol) { + unsigned int roundedTol = static_cast(llround(tol)); + return std::llabs(x - y) <= roundedTol; +} +template <> +DLLExport bool withinRelativeDifference(unsigned int const x, unsigned int const y, + double const tol) { + unsigned int roundedTol = static_cast(llround(tol)); + return std::llabs(x - y) <= roundedTol; +} + } // namespace Mantid::Kernel diff --git a/Framework/Kernel/src/Statistics.cpp b/Framework/Kernel/src/Statistics.cpp index 77e136158f5e..d638cba1beaa 100644 --- a/Framework/Kernel/src/Statistics.cpp +++ b/Framework/Kernel/src/Statistics.cpp @@ -92,6 +92,7 @@ template std::vector getZscore(const vector &data) } for (auto it = data.cbegin(); it != data.cend(); ++it) { auto tmp = static_cast(*it); + // unclear why Zscore is non-negative, was first implemented in #5316 Zscore.emplace_back(fabs((stats.mean - tmp) / stats.standard_deviation)); } return Zscore; diff --git a/Framework/Kernel/src/V3D.cpp b/Framework/Kernel/src/V3D.cpp index 64b0387018ef..d471546dc2ba 100644 --- a/Framework/Kernel/src/V3D.cpp +++ b/Framework/Kernel/src/V3D.cpp @@ -241,6 +241,7 @@ bool V3D::nullVector(const double tolerance) const noexcept { } bool V3D::unitVector(const double tolerance) const noexcept { + // NOTE could be made more efficient using norm2() const auto l = norm(); return std::abs(l - 1.) < tolerance; } @@ -500,12 +501,12 @@ V3D V3D::directionAngles(bool inDegrees) const { */ int V3D::maxCoeff() { int MaxOrder = 0; - if (abs(static_cast(m_pt[0])) > MaxOrder) - MaxOrder = abs(static_cast(m_pt[0])); - if (abs(static_cast(m_pt[1])) > MaxOrder) - MaxOrder = abs(static_cast(m_pt[1])); - if (abs(static_cast(m_pt[2])) > MaxOrder) - MaxOrder = abs(static_cast(m_pt[2])); + if (std::abs(static_cast(m_pt[0])) > MaxOrder) + MaxOrder = std::abs(static_cast(m_pt[0])); + if (std::abs(static_cast(m_pt[1])) > MaxOrder) + MaxOrder = std::abs(static_cast(m_pt[1])); + if (std::abs(static_cast(m_pt[2])) > MaxOrder) + MaxOrder = std::abs(static_cast(m_pt[2])); return MaxOrder; } @@ -513,15 +514,16 @@ int V3D::maxCoeff() { Calculates the absolute value. @return The absolute value */ -V3D V3D::absoluteValue() const { return V3D(fabs(m_pt[0]), fabs(m_pt[1]), fabs(m_pt[2])); } +V3D V3D::absoluteValue() const { return V3D(std::abs(m_pt[0]), std::abs(m_pt[1]), std::abs(m_pt[2])); } /** Calculates the error of the HKL to compare with tolerance @return The error */ double V3D::hklError() const { - return fabs(m_pt[0] - std::round(m_pt[0])) + fabs(m_pt[1] - std::round(m_pt[1])) + - fabs(m_pt[2] - std::round(m_pt[2])); + return std::abs(m_pt[0] - std::round(m_pt[0])) + // + std::abs(m_pt[1] - std::round(m_pt[1])) + // + std::abs(m_pt[2] - std::round(m_pt[2])); } } // namespace Mantid::Kernel diff --git a/Framework/Kernel/test/FloatingPointComparisonTest.h b/Framework/Kernel/test/FloatingPointComparisonTest.h index 9fff197fb606..39beaf690c72 100644 --- a/Framework/Kernel/test/FloatingPointComparisonTest.h +++ b/Framework/Kernel/test/FloatingPointComparisonTest.h @@ -16,15 +16,62 @@ class FloatingPointComparisonTest : public CxxTest::TestSuite { void test_Same_Value_Compare_Equal() { TS_ASSERT(Mantid::Kernel::equals(2.5, 2.5)); } void test_Difference_By_Machine_Eps_Compare_Equal() { - TS_ASSERT(Mantid::Kernel::equals(2.5, 2.5 + std::numeric_limits::epsilon())); + double const a = 0x1.4p1; // i.e. 2.5 + // increase by the machine precision + double const diff = std::ldexp(std::numeric_limits::epsilon(), std::ilogb(a)); + TS_ASSERT_DIFFERS(a, a + diff); + TS_ASSERT(Mantid::Kernel::equals(a, a + diff)); } void test_Difference_By_Machine_Eps_Plus_Small_Does_Not_Compare_Equal() { - TS_ASSERT_EQUALS(Mantid::Kernel::equals(2.5, 2.5 + 1.1 * std::numeric_limits::epsilon()), false); + double const a = 0x1.4p1; // i.e. 2.5 + // as above, but increase by twice the machine precision + double const diff = std::ldexp(std::numeric_limits::epsilon(), std::ilogb(a) + 1); + TS_ASSERT_DIFFERS(a, a + diff); + TS_ASSERT(!Mantid::Kernel::equals(a, a + diff)); + } + + void test_Similar_Small_Numbers_Compare_Equal() { + double const a = 0x1p-100; + // increase by the machine precision + double const diff = std::ldexp(std::numeric_limits::epsilon(), std::ilogb(a)); + TS_ASSERT_DIFFERS(a, a + diff); + TS_ASSERT(Mantid::Kernel::equals(a, a + diff)); + } + + void test_Different_Small_Numbers_Do_Not_Compare_Equal() { + // two small but machine-distinguishable numbers + double const a = 0x1.0p-100; // 1.0 * 2^{-100} + double const b = 0x1.8p-100; // 1.5 * 2^{-100} + double const diff = std::abs(a - b); + // the difference is less than machine epsilon (when scaled to 1) + TS_ASSERT_LESS_THAN(diff, std::numeric_limits::epsilon()); + // ne'ertheless, the numbers compare different + TS_ASSERT(!Mantid::Kernel::equals(a, b)) } void test_Same_Large_Numbers_Compare_Equal() { TS_ASSERT(Mantid::Kernel::equals(DBL_MAX, DBL_MAX)); } + void test_Similar_Large_Numbers_Compare_Equal() { + double const a = DBL_MAX / 2; + double const diff = std::ldexp(std::numeric_limits::epsilon(), std::ilogb(a)); + // the difference is a sizeable number and not by itself insignificant + TS_ASSERT_LESS_THAN(0x1.p50, diff); + // the numbers are technicaly different + TS_ASSERT_DIFFERS(a, a + diff); + // but they compare different + TS_ASSERT(Mantid::Kernel::equals(a, a + diff)); + } + + void test_Different_Large_Numbers_Do_Not_Compare_Equal() { + double const a = DBL_MAX / 2; + // as above, but increase by twice the machine precision + double const diff = std::ldexp(std::numeric_limits::epsilon(), std::ilogb(a) + 1); + TS_ASSERT_LESS_THAN(0x1.p50, diff); + TS_ASSERT_DIFFERS(a, a + diff); + TS_ASSERT(Mantid::Kernel::equals(a, a + diff)); + } + void test_Numbers_Outside_Custom_Tolerance_Are_Not_Equal() { const double tol(1e-08); TS_ASSERT_EQUALS(Mantid::Kernel::equals(0.1, 1.0001 * tol), false); diff --git a/buildconfig/CMake/CppCheck_Suppressions.txt.in b/buildconfig/CMake/CppCheck_Suppressions.txt.in index e6310c63f998..b906b1bc520d 100644 --- a/buildconfig/CMake/CppCheck_Suppressions.txt.in +++ b/buildconfig/CMake/CppCheck_Suppressions.txt.in @@ -170,10 +170,6 @@ unreadVariable:${CMAKE_SOURCE_DIR}/Framework/Algorithms/src/CalculateCarpenterSa missingOverride:${CMAKE_SOURCE_DIR}/Framework/Algorithms/inc/MantidAlgorithms/SpectrumAlgorithm.h:97 returnByReference:${CMAKE_SOURCE_DIR}/Framework/API/inc/MantidAPI/FunctionDomain1D.h:89 returnByReference:${CMAKE_SOURCE_DIR}/Framework/API/inc/MantidAPI/FunctionValues.h:86 -knownConditionTrueFalse:${CMAKE_SOURCE_DIR}/Framework/DataObjects/inc/MantidDataObjects/TableColumn.h:67 -knownConditionTrueFalse:${CMAKE_SOURCE_DIR}/Framework/DataObjects/inc/MantidDataObjects/TableColumn.h:70 -knownConditionTrueFalse:${CMAKE_SOURCE_DIR}/Framework/DataObjects/inc/MantidDataObjects/TableColumn.h:81 -knownConditionTrueFalse:${CMAKE_SOURCE_DIR}/Framework/DataObjects/inc/MantidDataObjects/TableColumn.h:84 constParameterReference:${CMAKE_SOURCE_DIR}/Framework/Algorithms/src/CalculatePolynomialBackground.cpp:174 missingOverride:${CMAKE_SOURCE_DIR}/Framework/Algorithms/inc/MantidAlgorithms/CalculatePlaczekSelfScattering.h:22 constVariableReference:${CMAKE_SOURCE_DIR}/Framework/Algorithms/src/ChangeTimeZero.cpp:154 @@ -202,7 +198,6 @@ constVariableReference:${CMAKE_SOURCE_DIR}/Framework/Algorithms/src/CreatePSDBle constVariablePointer:${CMAKE_SOURCE_DIR}/Framework/Algorithms/src/CreateLogPropertyTable.cpp:126 variableScope:${CMAKE_SOURCE_DIR}/Framework/Algorithms/src/CreateGroupingWorkspace.cpp:684 constVariablePointer:${CMAKE_SOURCE_DIR}/Framework/API/inc/MantidAPI/EnabledWhenWorkspaceIsType.h:50 -constVariableReference:${CMAKE_SOURCE_DIR}/Framework/Algorithms/src/DetectorEfficiencyCor.cpp:126 constVariablePointer:${CMAKE_SOURCE_DIR}/Framework/Algorithms/src/DetectorEfficiencyCorUser.cpp:176 constParameterReference:${CMAKE_SOURCE_DIR}/Framework/Algorithms/src/DirectILLTubeBackground.cpp:348 constParameterReference:${CMAKE_SOURCE_DIR}/Framework/Algorithms/src/DirectILLTubeBackground.cpp:364