Fix relative comparison in CompareWorkspaces - ornl-next

Framework/Algorithms/inc/MantidAlgorithms/CompareWorkspaces.h

@@ -78,6 +78,9 @@ 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);
+
 private:
   /// Initialise algorithm
   void init() override;
@@ -117,8 +120,6 @@ class MANTID_ALGORITHMS_DLL CompareWorkspaces final : public API::Algorithm {
   /// Records a mismatch in the Messages workspace and sets Result to false
   void recordMismatch(const std::string &msg, std::string ws1 = "", std::string ws2 = "");
 
-  bool relErr(double x1, double x2, double errorVal) const;
-
   /// Result of comparison (true if equal, false otherwise)
   bool m_result{false};
 

Framework/Algorithms/src/CompareWorkspaces.cpp

@@ -78,13 +78,13 @@ int compareEventLists(Kernel::Logger &logger, const EventList &el1, const EventL
       diffpulse = true;
       ++numdiffpulse;
     }
-    if (fabs(e1.tof() - e2.tof()) > tolTof) {
+    if (std::fabs(e1.tof() - e2.tof()) > tolTof) {
       difftof = true;
       ++numdifftof;
     }
     if (diffpulse && difftof)
       ++numdiffboth;
-    if (fabs(e1.weight() - e2.weight()) > tolWeight) {
+    if (std::fabs(e1.weight() - e2.weight()) > tolWeight) {
       diffweight = true;
       ++numdiffweight;
     }
@@ -265,7 +265,7 @@ void CompareWorkspaces::processGroups(const std::shared_ptr<const API::Workspace
     checker->setPropertyValue("Workspace1", namesOne[i]);
     checker->setPropertyValue("Workspace2", namesTwo[i]);
     for (size_t j = 0; j < numNonDefault; ++j) {
-      Property *p = nonDefaultProps[j];
+      Property const *p = nonDefaultProps[j];
       checker->setPropertyValue(p->name(), p->value());
     }
     checker->execute();
@@ -592,15 +592,26 @@ bool CompareWorkspaces::checkData(const API::MatrixWorkspace_const_sptr &ws1,
   // Cache a few things for later use
   const size_t numHists = ws1->getNumberHistograms();
   bool raggedWorkspace{false};
-  size_t numBins;
+  size_t numBins(0UL);
   try {
     numBins = ws1->blocksize();
   } catch (std::length_error &) {
     raggedWorkspace = true;
   }
   const bool histogram = ws1->isHistogramData();
   const bool checkAllData = getProperty("CheckAllData");
-  const bool RelErr = getProperty("ToleranceRelErr");
+  const bool isRelErr = getProperty("ToleranceRelErr");
+  const double tolerance = getProperty("Tolerance");
+  std::function<bool(double const, double const)> compare;
+  if (isRelErr) {
+    compare = [tolerance](double const x1, double const x2) -> bool {
+      return CompareWorkspaces::withinRelativeTolerance(x1, x2, tolerance);
+    };
+  } else {
+    compare = [tolerance](double const x1, double const x2) -> bool {
+      return CompareWorkspaces::withinAbsoluteTolerance(x1, x2, tolerance);
+    };
+  }
 
   // First check that the workspace are the same size
   if (numHists != ws2->getNumberHistograms() ||
@@ -615,8 +626,8 @@ bool CompareWorkspaces::checkData(const API::MatrixWorkspace_const_sptr &ws1,
     return false;
   }
 
-  const double tolerance = getProperty("Tolerance");
   bool resultBool = true;
+  bool logDebug = g_log.is(Logger::Priority::PRIO_DEBUG);
 
   // Now check the data itself
   PARALLEL_FOR_IF(m_parallelComparison && ws1->threadSafe() && ws2->threadSafe())
@@ -642,29 +653,26 @@ bool CompareWorkspaces::checkData(const API::MatrixWorkspace_const_sptr &ws1,
       } else {
 
         for (int j = 0; j < static_cast<int>(Y1.size()); ++j) {
-          bool err;
-          if (RelErr) {
-            err =
-                (relErr(X1[j], X2[j], tolerance) || relErr(Y1[j], Y2[j], tolerance) || relErr(E1[j], E2[j], tolerance));
-          } else
-            err = (std::fabs(X1[j] - X2[j]) > tolerance || std::fabs(Y1[j] - Y2[j]) > tolerance ||
-                   std::fabs(E1[j] - E2[j]) > tolerance);
-
+          bool err = (!compare(X1[j], X2[j]) || !compare(Y1[j], Y2[j]) || !compare(E1[j], E2[j]));
           if (err) {
-            g_log.debug() << "Data mismatch at cell (hist#,bin#): (" << i << "," << j << ")\n";
-            g_log.debug() << " Dataset #1 (X,Y,E) = (" << X1[j] << "," << Y1[j] << "," << E1[j] << ")\n";
-            g_log.debug() << " Dataset #2 (X,Y,E) = (" << X2[j] << "," << Y2[j] << "," << E2[j] << ")\n";
-            g_log.debug() << " Difference (X,Y,E) = (" << std::fabs(X1[j] - X2[j]) << "," << std::fabs(Y1[j] - Y2[j])
-                          << "," << std::fabs(E1[j] - E2[j]) << ")\n";
+            if (logDebug) {
+              g_log.debug() << "Data mismatch at cell (hist#,bin#): (" << i << "," << j << ")\n";
+              g_log.debug() << " Dataset #1 (X,Y,E) = (" << X1[j] << "," << Y1[j] << "," << E1[j] << ")\n";
+              g_log.debug() << " Dataset #2 (X,Y,E) = (" << X2[j] << "," << Y2[j] << "," << E2[j] << ")\n";
+              g_log.debug() << " Difference (X,Y,E) = (" << std::fabs(X1[j] - X2[j]) << "," << std::fabs(Y1[j] - Y2[j])
+                            << "," << std::fabs(E1[j] - E2[j]) << ")\n";
+            }
             PARALLEL_CRITICAL(resultBool)
             resultBool = false;
           }
         }
 
         // Extra one for histogram data
-        if (histogram && std::fabs(X1.back() - X2.back()) > tolerance) {
-          g_log.debug() << " Data ranges mismatch for spectra N: (" << i << ")\n";
-          g_log.debug() << " Last bin ranges (X1_end vs X2_end) = (" << X1.back() << "," << X2.back() << ")\n";
+        if (histogram && !compare(X1.back(), X2.back())) {
+          if (logDebug) {
+            g_log.debug() << " Data ranges mismatch for spectra N: (" << i << ")\n";
+            g_log.debug() << " Last bin ranges (X1_end vs X2_end) = (" << X1.back() << "," << X2.back() << ")\n";
+          }
           PARALLEL_CRITICAL(resultBool)
           resultBool = false;
         }
@@ -676,7 +684,7 @@ bool CompareWorkspaces::checkData(const API::MatrixWorkspace_const_sptr &ws1,
 
   if (!resultBool)
     recordMismatch("Data mismatch");
-  // If all is well, return true
+  // return result
   return resultBool;
 }
 
@@ -944,7 +952,7 @@ bool CompareWorkspaces::checkRunProperties(const API::Run &run1, const API::Run
     return false;
   } else {
     // Sort logs by name before one-by-one comparison
-    auto compareNames = [](Kernel::Property *p1, Kernel::Property *p2) { return p1->name() < p2->name(); };
+    auto compareNames = [](Kernel::Property const *p1, Kernel::Property const *p2) { return p1->name() < p2->name(); };
     std::sort(ws1logs.begin(), ws1logs.end(), compareNames);
     std::sort(ws2logs.begin(), ws2logs.end(), compareNames);
     for (size_t i = 0; i < ws1logs.size(); ++i) {
@@ -1097,17 +1105,13 @@ void CompareWorkspaces::doPeaksComparison(PeaksWorkspace_sptr tws1, PeaksWorkspa
       }
       bool mismatch = false;
       if (isRelErr) {
-        if (relErr(s1, s2, tolerance)) {
-          mismatch = true;
-        }
-      } else if (std::fabs(s1 - s2) > tolerance) {
-        mismatch = true;
-      } else if (std::fabs(v1[0] - v2[0]) > tolerance) {
-        mismatch = true;
-      } else if (std::fabs(v1[1] - v2[1]) > tolerance) {
-        mismatch = true;
-      } else if (std::fabs(v1[2] - v2[2]) > tolerance) {
-        mismatch = true;
+        mismatch = !withinRelativeTolerance(s1, s2, tolerance);
+        // Q: whould we not also compare the vectors?
+      } else {
+        mismatch = !withinAbsoluteTolerance(s1, s2, tolerance) ||       //
+                   !withinAbsoluteTolerance(v1[0], v2[0], tolerance) || //
+                   !withinAbsoluteTolerance(v1[1], v2[1], tolerance) || //
+                   !withinAbsoluteTolerance(v1[2], v2[2], tolerance);   //
       }
       if (mismatch) {
         g_log.notice(name);
@@ -1213,11 +1217,12 @@ void CompareWorkspaces::doLeanElasticPeaksComparison(const LeanElasticPeaksWorks
         g_log.information() << "Column " << name << " is not compared\n";
       }
       bool mismatch = false;
+      // Q: why does it not perform the user-specified operation for QLab and QSample?
       if (isRelErr && name != "QLab" && name != "QSample") {
-        if (relErr(s1, s2, tolerance)) {
+        if (!withinRelativeTolerance(s1, s2, tolerance)) {
           mismatch = true;
         }
-      } else if (std::fabs(s1 - s2) > tolerance) {
+      } else if (!withinAbsoluteTolerance(s1, s2, tolerance)) {
         mismatch = true;
       }
       if (mismatch) {
@@ -1275,19 +1280,13 @@ void CompareWorkspaces::doTableComparison(const API::ITableWorkspace_const_sptr
     const auto c2 = tws2->getColumn(i);
 
     if (isRelErr) {
-      if (!c1->equalsRelErr(*c2, tolerance)) {
-        mismatch = true;
-      }
+      mismatch = !c1->equalsRelErr(*c2, tolerance);
     } else {
-
-      if (!c1->equals(*c2, tolerance)) {
-        mismatch = true;
-      }
+      mismatch = !c1->equals(*c2, tolerance);
     }
     if (mismatch) {
       g_log.debug() << "Table data mismatch at column " << i << "\n";
       recordMismatch("Table data mismatch");
-      mismatch = false;
       if (!checkAllData) {
         return;
       }
@@ -1340,27 +1339,41 @@ void CompareWorkspaces::recordMismatch(const std::string &msg, std::string ws1,
   m_result = false;
 }
 
+//------------------------------------------------------------------------------------------------
+/** Function which calculates absolute error between two values and analyses if
+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
+
+@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::fabs(x1 - x2) > atol);
+}
+
 //------------------------------------------------------------------------------------------------
 /** Function which calculates relative error between two values and analyses if
-this error is within the limits
-* requested. When the absolute value of the difference is smaller then the value
-of the error requested,
-* absolute error is used instead of relative error.
+this error is within the limits requested.
 
-@param x1       -- first value to check difference
-@param x2       -- second value to check difference
-@param errorVal -- the value of the error, to check against. Should  be large
-then 0
+@param x1    -- first value to check difference
+@param x2    -- second value to check difference
+@param rtol  -- the tolerance of the comparison. Must be nonnegative
 
-@returns true if error or false if the value is within the limits requested
+@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::relErr(double x1, double x2, double errorVal) const {
-  double num = std::fabs(x1 - x2);
-  // how to treat x1<0 and x2 > 0 ?  probably this way
-  double den = 0.5 * (std::fabs(x1) + std::fabs(x2));
-  if (den < errorVal)
-    return (num > errorVal);
-
-  return (num / den > errorVal);
+bool CompareWorkspaces::withinRelativeTolerance(double const x1, double const x2, double const rtol) {
+  // calculate difference
+  double const num = std::fabs(x1 - x2);
+  // return early if the values are equal
+  if (num == 0.0)
+    return true;
+  // compare the difference to the midpoint value -- could lead to issues for negative values
+  double const den = 0.5 * std::fabs(x1 + x2);
+  // NOTE !(num > rtol*den) is not the same as (num <= rtol*den)
+  return !(num > (rtol * den));
 }
 } // namespace Mantid::Algorithms

Framework/Algorithms/test/CreateUserDefinedBackgroundTest.h

@@ -23,6 +23,9 @@ using Mantid::API::ITableWorkspace_sptr;
 using Mantid::API::MatrixWorkspace_sptr;
 
 namespace {
+// flag for absolute or relative comparison
+enum Comparison : bool { CMP_RELATIVE = true, CMP_ABSOLUTE = false };
+
 // Gaussian
 double gauss(double x, double height, double centre, double fwhm) {
   const double factor = 2.0 * sqrt(2.0 * log(2.0));
@@ -108,7 +111,7 @@ class CreateUserDefinedBackgroundTest : public CxxTest::TestSuite {
 
     // The expected result
     const auto expected = createExpectedResults(true);
-    TS_ASSERT(workspacesEqual(expected, outputWS, 0.1, true));
+    TS_ASSERT(workspacesEqual(expected, outputWS, 0.105, Comparison::CMP_RELATIVE));
   }
 
   void test_exec_HistoWS_NormalisePlotsOn() {
@@ -134,7 +137,7 @@ class CreateUserDefinedBackgroundTest : public CxxTest::TestSuite {
 
     // The expected result
     const auto expected = createExpectedResults(true, true);
-    TS_ASSERT(workspacesEqual(expected, outputWS, 5e-2));
+    TS_ASSERT(workspacesEqual(expected, outputWS, 5e-2, Comparison::CMP_ABSOLUTE));
   }
 
   void test_exec_PointsWS_extend() {
@@ -244,7 +247,7 @@ class CreateUserDefinedBackgroundTest : public CxxTest::TestSuite {
     const auto expected = createExpectedResults(true, false);
     Mantid::API::WorkspaceHelpers::makeDistribution(expected);
 
-    TS_ASSERT(workspacesEqual(expected, outputWS, 0.1, true));
+    TS_ASSERT(workspacesEqual(expected, outputWS, 0.105, Comparison::CMP_RELATIVE));
   }
 
 private:
@@ -302,7 +305,7 @@ class CreateUserDefinedBackgroundTest : public CxxTest::TestSuite {
 
   /// Compare workspaces
   bool workspacesEqual(const MatrixWorkspace_sptr &lhs, const MatrixWorkspace_sptr &rhs, double tolerance,
-                       bool relativeError = false) {
+                       bool const relativeError = Comparison::CMP_ABSOLUTE) {
     auto alg = Mantid::API::AlgorithmFactory::Instance().create("CompareWorkspaces", 1);
     alg->setChild(true);
     alg->initialize();
@@ -336,7 +339,7 @@ class CreateUserDefinedBackgroundTest : public CxxTest::TestSuite {
 
     // The expected result
     const auto expected = createExpectedResults(false);
-    TS_ASSERT(workspacesEqual(expected, outputWS, 1e-4));
+    TS_ASSERT(workspacesEqual(expected, outputWS, 1e-4, Comparison::CMP_ABSOLUTE));
   }
 
   /// Cached string for option

Testing/SystemTests/tests/framework/DirectILLAutoProcessTest.py

@@ -29,7 +29,7 @@ def cleanup(self):
 
     def validate(self):
         self.tolerance = 1e-2
-        self.tolerance_is_rel_err = True
+        self.tolerance_is_rel_err = False
         self.disableChecking = ["Instrument", "SpectraMap"]
         return ["He3C60", "ILL_PANTHER_Powder_Auto.nxs"]
 

Testing/SystemTests/tests/framework/ILLDirectGeometryReductionTest.py

@@ -164,7 +164,7 @@ def runTest(self):
 
     def validate(self):
         self.tolerance = 1e-2
-        self.tolerance_is_rel_err = True
+        self.tolerance_is_rel_err = False
         self.disableChecking = ["Instrument", "Sample"]
         return ["cropped", "ILL_IN6_SofQW.nxs"]
 

Testing/SystemTests/tests/framework/PaalmanPingsMonteCarloAbsorptionTest.py

@@ -91,7 +91,7 @@ def cleanup(self):
 
     def validate(self):
         self.tolerance = 1e-3
-        self.tolerance_is_rel_err = True
+        self.tolerance_is_rel_err = False
         return ["annulus_corr", "irs_PP_MC_annulus.nxs"]
 
     def runTest(self):

buildconfig/CMake/CppCheck_Suppressions.txt.in

@@ -204,8 +204,6 @@ constVariableReference:${CMAKE_SOURCE_DIR}/Framework/Algorithms/src/CreatePSDBle
 constVariableReference:${CMAKE_SOURCE_DIR}/Framework/Algorithms/src/CreatePSDBleedMask.cpp:214
 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/Algorithms/src/CompareWorkspaces.cpp:268
-constParameterPointer:${CMAKE_SOURCE_DIR}/Framework/Algorithms/src/CompareWorkspaces.cpp:947
 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
@@ -1676,3 +1674,4 @@ returnByReference:${CMAKE_SOURCE_DIR}/qt/scientific_interfaces/Indirect/Reductio
 returnByReference:${CMAKE_SOURCE_DIR}/qt/scientific_interfaces/Indirect/Reduction/ISISEnergyTransferData.h:193
 returnByReference:${CMAKE_SOURCE_DIR}/qt/scientific_interfaces/Indirect/Reduction/ISISEnergyTransferData.h:199
 returnByReference:${CMAKE_SOURCE_DIR}/qt/scientific_interfaces/Indirect/Reduction/ISISEnergyTransferData.h:218
+unreadVariable:${CMAKE_SOURCE_DIR}/Framework/Algorithms/src/CompareWorkspaces.cpp:1277

docs/source/release/v6.11.0/Indirect/Algorithms/Bugfixes/37889.rst

@@ -0,0 +1 @@
+- fix :ref:`CompareWorkspaces <algm-CompareWorkspaces-v1>` for relative differences of small values.