replace certain parts of algorithms with floating point ops

Framework/API/src/MatrixWorkspace.cpp

@@ -1951,7 +1951,7 @@ std::pair<size_t, double> MatrixWorkspace::getXIndex(size_t i, double x, bool is
       auto index = static_cast<size_t>(std::distance(X.begin(), ix));
       if (isLeft)
         --index;
-      return std::make_pair(index, fabs((X[index] - x) / (*ix - *(ix - 1))));
+      return std::make_pair(index, std::abs((X[index] - x) / (*ix - *(ix - 1))));
     }
   }
   // I don't think we can ever get here

Framework/API/src/NumericAxis.cpp

@@ -9,6 +9,7 @@
 //----------------------------------------------------------------------
 #include "MantidAPI/NumericAxis.h"
 #include "MantidKernel/Exception.h"
+#include "MantidKernel/FloatingPointComparison.h"
 #include "MantidKernel/VectorHelper.h"
 
 #include <boost/format.hpp>
@@ -22,12 +23,17 @@ Mantid::Kernel::Logger g_log("NumericAxis");
 class EqualWithinTolerance {
 public:
   explicit EqualWithinTolerance(double tolerance) : m_tolerance(tolerance){};
+  /**
+   * This handles NaNs and infs differently than the FloatingPointComparison operations.
+   * If this is not necessary, then this entire class may be replaced with
+   * Mantid::Kernel::withinAbsoluteDifference
+   */
   bool operator()(double a, double b) {
     if (std::isnan(a) && std::isnan(b))
       return true;
     if (std::isinf(a) && std::isinf(b))
       return true;
-    return std::abs(a - b) <= m_tolerance;
+    return Mantid::Kernel::withinAbsoluteDifference(a, b, m_tolerance);
   }
 
 private:

Framework/API/src/WorkspaceOpOverloads.cpp

@@ -13,6 +13,7 @@
 #include "MantidAPI/MatrixWorkspace.h"
 #include "MantidAPI/WorkspaceFactory.h"
 #include "MantidAPI/WorkspaceGroup.h"
+#include "MantidKernel/FloatingPointComparison.h"
 #include "MantidKernel/Property.h"
 
 #include <numeric>
@@ -376,10 +377,10 @@ bool WorkspaceHelpers::matchingBins(const MatrixWorkspace &ws1, const MatrixWork
   const double secondWS = std::accumulate(ws2.x(0).begin(), ws2.x(0).end(), 0.);
   if (std::abs(firstWS) < 1.0E-7 && std::abs(secondWS) < 1.0E-7) {
     for (size_t i = 0; i < ws1.x(0).size(); i++) {
-      if (std::abs(ws1.x(0)[i] - ws2.x(0)[i]) > 1.0E-7)
+      if (!Kernel::withinAbsoluteDifference(ws1.x(0)[i], ws2.x(0)[i], 1.0E-7))
         return false;
     }
-  } else if (std::abs(firstWS - secondWS) / std::max<double>(std::abs(firstWS), std::abs(secondWS)) > 1.0E-7)
+  } else if (!Kernel::withinRelativeDifference(firstWS, secondWS, 1.0E-7))
     return false;
 
   // If we were only asked to check the first spectrum, return now
@@ -409,11 +410,10 @@ bool WorkspaceHelpers::matchingBins(const MatrixWorkspace &ws1, const MatrixWork
     const double secondWSLoop = std::accumulate(ws2.x(i).begin(), ws2.x(i).end(), 0.);
     if (std::abs(firstWSLoop) < 1.0E-7 && std::abs(secondWSLoop) < 1.0E-7) {
       for (size_t j = 0; j < ws1.x(i).size(); j++) {
-        if (std::abs(ws1.x(i)[j] - ws2.x(i)[j]) > 1.0E-7)
+        if (!Kernel::withinAbsoluteDifference(ws1.x(i)[j], ws2.x(i)[j], 1.0E-7))
           return false;
       }
-    } else if (std::abs(firstWSLoop - secondWSLoop) / std::max<double>(std::abs(firstWSLoop), std::abs(secondWSLoop)) >
-               1.0E-7)
+    } else if (!Kernel::withinRelativeDifference(firstWSLoop, secondWSLoop, 1.0E-7))
       return false;
   }
 

Framework/Algorithms/src/CalculateCarpenterSampleCorrection.cpp

@@ -15,6 +15,7 @@
 #include "MantidDataObjects/WorkspaceCreation.h"
 #include "MantidGeometry/Instrument.h"
 #include "MantidKernel/CompositeValidator.h"
+#include "MantidKernel/FloatingPointComparison.h"
 #include "MantidKernel/Material.h"
 
 #include <stdexcept>
@@ -128,12 +129,11 @@ void CalculateCarpenterSampleCorrection::exec() {
   const Material &sampleMaterial = inputWksp->sample().getMaterial();
   if (sampleMaterial.totalScatterXSection() != 0.0) {
     g_log.information() << "Using material \"" << sampleMaterial.name() << "\" from workspace\n";
-    if (std::abs(coeff1 - 2.8) < std::numeric_limits<double>::epsilon())
+    if (Kernel::equals(coeff1, 2.8))
       coeff1 = sampleMaterial.absorbXSection(LAMBDA_REF) / LAMBDA_REF;
-    if ((std::abs(coeff2 - 0.0721) < std::numeric_limits<double>::epsilon()) &&
-        (!isEmpty(sampleMaterial.numberDensity())))
+    if (Kernel::equals(coeff2, 0.0721) && !isEmpty(sampleMaterial.numberDensity()))
       coeff2 = sampleMaterial.numberDensity();
-    if (std::abs(coeff3 - 5.1) < std::numeric_limits<double>::epsilon())
+    if (Kernel::equals(coeff3, 5.1))
       coeff3 = sampleMaterial.totalScatterXSection();
   } else // Save input in Sample with wrong atomic number and name
   {

Framework/Algorithms/src/DetectorEfficiencyCor.cpp

@@ -19,6 +19,7 @@
 #include "MantidKernel/BoundedValidator.h"
 #include "MantidKernel/CompositeValidator.h"
 #include "MantidKernel/Exception.h"
+#include "MantidKernel/FloatingPointComparison.h"
 #include "MantidKernel/PhysicalConstants.h"
 #include "MantidTypes/SpectrumDefinition.h"
 
@@ -301,7 +302,7 @@ void DetectorEfficiencyCor::getDetectorGeometry(const Geometry::IDetector &det,
   if (it == m_shapeCache.end()) {
     double xDist = distToSurface(V3D(DIST_TO_UNIVERSE_EDGE, 0, 0), shape_sptr.get());
     double zDist = distToSurface(V3D(0, 0, DIST_TO_UNIVERSE_EDGE), shape_sptr.get());
-    if (std::abs(zDist - xDist) < 1e-8) {
+    if (withinAbsoluteDifference(zDist, xDist, 1e-6)) {
       detRadius = zDist / 2.0;
       detAxis = V3D(0, 1, 0);
       // assume radi in z and x and the axis is in the y
@@ -311,7 +312,7 @@ void DetectorEfficiencyCor::getDetectorGeometry(const Geometry::IDetector &det,
       return;
     }
     double yDist = distToSurface(V3D(0, DIST_TO_UNIVERSE_EDGE, 0), shape_sptr.get());
-    if (std::abs(yDist - zDist) < 1e-8) {
+    if (withinAbsoluteDifference(yDist, zDist, 1e-8)) {
       detRadius = yDist / 2.0;
       detAxis = V3D(1, 0, 0);
       // assume that y and z are radi of the cylinder's circular cross-section
@@ -322,7 +323,7 @@ void DetectorEfficiencyCor::getDetectorGeometry(const Geometry::IDetector &det,
       return;
     }
 
-    if (std::abs(xDist - yDist) < 1e-8) {
+    if (withinAbsoluteDifference(xDist, yDist, 1e-8)) {
       detRadius = xDist / 2.0;
       detAxis = V3D(0, 0, 1);
       PARALLEL_CRITICAL(deteff_shapecachec) {

Framework/Algorithms/src/FindCenterOfMassPosition2.cpp

@@ -18,6 +18,7 @@
 #include "MantidKernel/ArrayProperty.h"
 #include "MantidKernel/BoundedValidator.h"
 #include "MantidKernel/CompositeValidator.h"
+#include "MantidKernel/FloatingPointComparison.h"
 #include "MantidKernel/PhysicalConstants.h"
 namespace Mantid::Algorithms {
 
@@ -29,11 +30,6 @@ using namespace API;
 using namespace Geometry;
 using namespace DataObjects;
 
-namespace { // anonymous namespace
-/// Equal within machine tolerance
-bool almost_equals(const double a, const double b) { return fabs(a - b) < std::numeric_limits<double>::min(); }
-} // anonymous namespace
-
 void FindCenterOfMassPosition2::init() {
   const auto wsValidator = std::make_shared<CompositeValidator>();
   const auto positiveDouble = std::make_shared<BoundedValidator<double>>();
@@ -113,7 +109,7 @@ void FindCenterOfMassPosition2::findCenterOfMass(const API::MatrixWorkspace_sptr
 
     // Check to see if we have the same result
     // as the previous iteration
-    if (almost_equals(distanceFromPrevious, distanceFromPreviousPrevious)) {
+    if (Kernel::equals(distanceFromPrevious, distanceFromPreviousPrevious)) {
       totalLocalMinima++;
     } else {
       totalLocalMinima = 0;

Framework/Algorithms/src/FitPeak.cpp

@@ -237,7 +237,7 @@ void FitOneSinglePeak::setupGuessedFWHM(double usrwidth, int minfwhm, int maxfwh
 void FitOneSinglePeak::setFitPeakCriteria(bool usepeakpostol, double peakpostol) {
   m_usePeakPositionTolerance = usepeakpostol;
   if (usepeakpostol) {
-    m_peakPositionTolerance = fabs(peakpostol);
+    m_peakPositionTolerance = std::abs(peakpostol);
     if (peakpostol < 1.0E-13)
       g_log.warning("Peak position tolerance is very tight. ");
   }
@@ -931,7 +931,7 @@ void FitOneSinglePeak::processNStoreFitResult(double rwp, bool storebkgd) {
     double f_centre = m_peakFunc->centre();
     if (m_usePeakPositionTolerance) {
       // Peak position criteria is on position tolerance
-      if (fabs(f_centre - m_userPeakCentre) > m_peakPositionTolerance) {
+      if (!Kernel::withinAbsoluteDifference(f_centre, m_userPeakCentre, m_peakPositionTolerance)) {
         rwp = DBL_MAX;
         failreason = "Peak centre out of tolerance. ";
         fitsuccess = false;

Framework/Algorithms/src/He3TubeEfficiency.cpp

@@ -20,6 +20,7 @@
 #include "MantidKernel/ArrayBoundedValidator.h"
 #include "MantidKernel/ArrayProperty.h"
 #include "MantidKernel/CompositeValidator.h"
+#include "MantidKernel/FloatingPointComparison.h"
 
 #include <cmath>
 #include <stdexcept>
@@ -216,7 +217,7 @@ double He3TubeEfficiency::calculateExponential(std::size_t spectraIndex, const G
   double sinTheta = std::sqrt(1.0 - cosTheta * cosTheta);
 
   const double straight_path = detDiameter - twiceTubeThickness;
-  if (std::fabs(straight_path - 0.0) < TOL) {
+  if (Kernel::withinAbsoluteDifference(straight_path, 0.0, TOL)) {
     throw std::out_of_range("Twice tube thickness cannot be greater than "
                             "or equal to the tube diameter");
   }
@@ -245,15 +246,15 @@ void He3TubeEfficiency::getDetectorGeometry(const Geometry::IDetector &det, doub
   if (it == m_shapeCache.end()) {
     double xDist = distToSurface(Kernel::V3D(DIST_TO_UNIVERSE_EDGE, 0, 0), shape_sptr.get());
     double zDist = distToSurface(Kernel::V3D(0, 0, DIST_TO_UNIVERSE_EDGE), shape_sptr.get());
-    if (std::abs(zDist - xDist) < 1e-8) {
+    if (Kernel::withinAbsoluteDifference(zDist, xDist, 1e-8)) {
       detRadius = zDist / 2.0;
       detAxis = Kernel::V3D(0, 1, 0);
       // assume radii in z and x and the axis is in the y
       PARALLEL_CRITICAL(deteff_shapecachea) { m_shapeCache.insert({shape_sptr.get(), {detRadius, detAxis}}); }
       return;
     }
     double yDist = distToSurface(Kernel::V3D(0, DIST_TO_UNIVERSE_EDGE, 0), shape_sptr.get());
-    if (std::abs(yDist - zDist) < 1e-8) {
+    if (Kernel::withinAbsoluteDifference(yDist, zDist, 1e-8)) {
       detRadius = yDist / 2.0;
       detAxis = Kernel::V3D(1, 0, 0);
       // assume that y and z are radii of the cylinder's circular cross-section
@@ -262,7 +263,7 @@ void He3TubeEfficiency::getDetectorGeometry(const Geometry::IDetector &det, doub
       return;
     }
 
-    if (std::abs(xDist - yDist) < 1e-8) {
+    if (Kernel::withinAbsoluteDifference(xDist, yDist, 1e-8)) {
       detRadius = xDist / 2.0;
       detAxis = Kernel::V3D(0, 0, 1);
       PARALLEL_CRITICAL(deteff_shapecachec) { m_shapeCache.insert({shape_sptr.get(), {detRadius, detAxis}}); }

Framework/Algorithms/src/Integration.cpp

@@ -18,6 +18,7 @@
 #include "MantidHistogramData/Histogram.h"
 #include "MantidKernel/ArrayProperty.h"
 #include "MantidKernel/BoundedValidator.h"
+#include "MantidKernel/FloatingPointComparison.h"
 #include "MantidKernel/VectorHelper.h"
 
 #include <cmath>
@@ -68,7 +69,7 @@ struct tolerant_less {
   bool operator()(const double &left, const double &right) const {
     // soft equal, if the diff left-right is below a numerical error
     // (uncertainty) threshold, we cannot say
-    return (left < right) && (std::abs(left - right) > 1 * std::numeric_limits<double>::epsilon());
+    return (left < right) && !Kernel::equals(left, right);
   }
 };
 

Framework/Algorithms/src/RealFFT.cpp

@@ -79,6 +79,8 @@ void RealFFT::exec() {
   double dx = (X.back() - X.front()) / static_cast<double>(X.size() - 1);
   if (!IgnoreXBins) {
     for (size_t i = 0; i < X.size() - 2; i++)
+      // note this cannot be replaced with Kernel::withinRelativeDifference,
+      // or fails to detect some errors
       if (std::abs(dx - X[i + 1] + X[i]) / dx > 1e-7)
         throw std::invalid_argument("X axis must be linear (all bins have same "
                                     "width). This can be ignored if "

Framework/Algorithms/src/RunCombinationHelpers/SampleLogsBehaviour.cpp

@@ -10,6 +10,7 @@
 #include "MantidAPI/Run.h"
 #include "MantidAlgorithms/RunCombinationHelpers/SampleLogsBehaviour.h"
 #include "MantidGeometry/Instrument.h"
+#include "MantidKernel/FloatingPointComparison.h"
 #include "MantidKernel/Property.h"
 #include "MantidKernel/StringTokenizer.h"
 #include "MantidKernel/Strings.h"
@@ -531,7 +532,7 @@ void SampleLogsBehaviour::checkErrorProperty(const MatrixWorkspace &addeeWS, Pro
 bool SampleLogsBehaviour::isWithinTolerance(const SampleLogBehaviour &behaviour, const double addeeWSNumericValue,
                                             const double outWSNumericValue) {
   if (behaviour.isNumeric && behaviour.tolerance > 0.0) {
-    return std::abs(addeeWSNumericValue - outWSNumericValue) < behaviour.tolerance;
+    return Kernel::withinAbsoluteDifference(addeeWSNumericValue, outWSNumericValue, behaviour.tolerance);
   }
 
   return false;

Framework/Algorithms/src/Stitch1D.cpp

@@ -16,6 +16,7 @@
 #include "MantidHistogramData/HistogramY.h"
 #include "MantidKernel/ArrayProperty.h"
 #include "MantidKernel/BoundedValidator.h"
+#include "MantidKernel/FloatingPointComparison.h"
 #include "MantidKernel/MultiThreaded.h"
 #include "MantidKernel/PropertyWithValue.h"
 #include "MantidKernel/RebinParamsValidator.h"
@@ -491,10 +492,10 @@ void Stitch1D::exec() {
     const double xMin = params.front();
     const double xMax = params.back();
 
-    if (std::abs(xMin - startOverlap) < 1E-6)
+    if (Kernel::withinAbsoluteDifference(xMin, startOverlap, 1E-6))
       startOverlap = xMin;
 
-    if (std::abs(xMax - endOverlap) < 1E-6)
+    if (Kernel::withinAbsoluteDifference(xMax, endOverlap, 1E-6))
       endOverlap = xMax;
 
     if (startOverlap < xMin) {

Framework/Algorithms/src/UnwrapMonitor.cpp

@@ -16,6 +16,7 @@
 #include "MantidGeometry/Instrument.h"
 #include "MantidKernel/BoundedValidator.h"
 #include "MantidKernel/CompositeValidator.h"
+#include "MantidKernel/FloatingPointComparison.h"
 #include "MantidKernel/PhysicalConstants.h"
 #include "MantidKernel/UnitFactory.h"
 
@@ -204,7 +205,7 @@ const std::vector<int> UnwrapMonitor::unwrapX(std::vector<double> &newX, const i
       const double wavelength = m_conversionConstant / velocity;
       newX.emplace_back(wavelength);
       // Remove the duplicate boundary bin
-      if (tof == m_Tmax && std::abs(wavelength - tempX_L.front()) < 1.0e-5)
+      if (tof == m_Tmax && Kernel::withinAbsoluteDifference(wavelength, tempX_L.front(), 1.0e-5))
         newX.pop_back();
       // Record the bins that fall in this range for copying over the data &
       // errors
@@ -223,7 +224,7 @@ const std::vector<int> UnwrapMonitor::unwrapX(std::vector<double> &newX, const i
 
   // Record the point at which the unwrapped sections are joined, first time
   // through only
-  Property *join = getProperty("JoinWavelength");
+  Property const *join = getProperty("JoinWavelength");
   if (join->isDefault()) {
     g_log.information() << "Joining wavelength: " << tempX_L.front() << " Angstrom\n";
     setProperty("JoinWavelength", tempX_L.front());

Framework/Crystal/src/FindSXPeaksHelper.cpp

@@ -8,6 +8,7 @@
 #include "MantidAPI/Progress.h"
 #include "MantidGeometry/Instrument/DetectorGroup.h"
 #include "MantidKernel/ConfigService.h"
+#include "MantidKernel/FloatingPointComparison.h"
 #include "MantidKernel/Logger.h"
 #include "MantidKernel/PhysicalConstants.h"
 #include "MantidKernel/Unit.h"
@@ -118,11 +119,11 @@ Object comparision
 @param tolerance : tolerance
 */
 bool SXPeak::compare(const SXPeak &rhs, double tolerance) const {
-  if (std::abs(m_tof / m_nPixels - rhs.m_tof / rhs.m_nPixels) > tolerance * m_tof / m_nPixels)
+  if (!Kernel::withinRelativeDifference(m_tof / m_nPixels, rhs.m_tof / rhs.m_nPixels, tolerance))
     return false;
-  if (std::abs(m_phi / m_nPixels - rhs.m_phi / rhs.m_nPixels) > tolerance * m_phi / m_nPixels)
+  if (!Kernel::withinRelativeDifference(m_phi / m_nPixels, rhs.m_phi / rhs.m_nPixels, tolerance))
     return false;
-  if (std::abs(m_twoTheta / m_nPixels - rhs.m_twoTheta / rhs.m_nPixels) > tolerance * m_twoTheta / m_nPixels)
+  if (!Kernel::withinRelativeDifference(m_twoTheta / m_nPixels, rhs.m_twoTheta / rhs.m_nPixels, tolerance))
     return false;
   return true;
 }

Framework/Crystal/src/PredictPeaks.cpp

@@ -25,6 +25,7 @@
 #include "MantidGeometry/Objects/InstrumentRayTracer.h"
 #include "MantidKernel/BoundedValidator.h"
 #include "MantidKernel/EnabledWhenProperty.h"
+#include "MantidKernel/FloatingPointComparison.h"
 #include "MantidKernel/ListValidator.h"
 
 #include <fstream>
@@ -201,7 +202,7 @@ void PredictPeaks::exec() {
     // Get all goniometer matrices
     DblMatrix lastGoniometerMatrix = Matrix<double>(3, 3, false);
     for (int i = 0; i < static_cast<int>(peaksWS->getNumberPeaks()); ++i) {
-      IPeak &p = peaksWS->getPeak(i);
+      IPeak const &p = peaksWS->getPeak(i);
       DblMatrix currentGoniometerMatrix = p.getGoniometerMatrix();
       if (!(currentGoniometerMatrix == lastGoniometerMatrix)) {
         gonioVec.emplace_back(currentGoniometerMatrix);
@@ -289,7 +290,7 @@ void PredictPeaks::exec() {
     m_detectorCacheSearch = std::make_unique<DetectorSearcher>(m_inst, m_pw->detectorInfo());
 
   if (!usingInstrument) {
-    for (auto &possibleHKL : possibleHKLs) {
+    for (auto const &possibleHKL : possibleHKLs) {
       calculateQAndAddToOutputLeanElastic(possibleHKL, ub);
     }
   } else if (getProperty("CalculateGoniometerForCW")) {
@@ -329,7 +330,7 @@ void PredictPeaks::exec() {
       if (!std::isfinite(angle) || angle < angleMin || angle > angleMax)
         continue;
 
-      if (std::abs(wavelength - lambda) < 0.01) {
+      if (Kernel::withinAbsoluteDifference(wavelength, lambda, 0.01)) {
         g_log.information() << "Found goniometer rotation to be in YZY convention [" << angles[0] << ", " << angles[1]
                             << ". " << angles[2] << "] degrees for Q sample = " << q_sample << "\n";
         calculateQAndAddToOutput(possibleHKL, orientedUB, goniometer.getR());
@@ -341,7 +342,7 @@ void PredictPeaks::exec() {
     logNumberOfPeaksFound(allowedPeakCount);
 
   } else {
-    for (auto &goniometerMatrix : gonioVec) {
+    for (auto const &goniometerMatrix : gonioVec) {
       // Final transformation matrix (HKL to Q in lab frame)
       DblMatrix orientedUB = goniometerMatrix * ub;
 
@@ -357,7 +358,7 @@ void PredictPeaks::exec() {
                            "no extended detector space has been defined\n";
       }
 
-      for (auto &possibleHKL : possibleHKLs) {
+      for (auto const &possibleHKL : possibleHKLs) {
         if (lambdaFilter.isAllowed(possibleHKL)) {
           calculateQAndAddToOutput(possibleHKL, orientedUB, goniometerMatrix);
           ++allowedPeakCount;
@@ -503,7 +504,7 @@ void PredictPeaks::fillPossibleHKLsUsingPeaksWorkspace(const IPeaksWorkspace_spt
   double peaks_q_convention_factor = qConventionFactor(peaksWorkspace->getConvention());
 
   for (int i = 0; i < static_cast<int>(peaksWorkspace->getNumberPeaks()); ++i) {
-    IPeak &p = peaksWorkspace->getPeak(i);
+    IPeak const &p = peaksWorkspace->getPeak(i);
     // Get HKL from that peak
     V3D hkl = p.getHKL() * peaks_q_convention_factor;