improves numerical test stability

This improves the stability of the tests and makes them independent of the utilized OS/ Compiler/ Architecture by:
- adding dedicated methods for FloatingPoint comparison using
  - UPL distance
  - Relative Epsilon
- adding custom GoogleTestMatchers using these new methods
- updating the DataStructures Equality methods

src/polyhedralGravity/model/GravityModelData.h

@@ -84,10 +84,10 @@ namespace polyhedralGravity {
          * @note Just used for testing purpose
          */
         bool operator==(const Distance &rhs) const {
-            return util::floatNear(l1, rhs.l1) &&
-                   util::floatNear(l2, rhs.l2) &&
-                   util::floatNear(s1, rhs.s1) &&
-                   util::floatNear(s2, rhs.s2);
+            return util::almostEqualRelative(l1, rhs.l1) &&
+                   util::almostEqualRelative(l2, rhs.l2) &&
+                   util::almostEqualRelative(s1, rhs.s1) &&
+                   util::almostEqualRelative(s2, rhs.s2);
         }
 
         /**
@@ -141,7 +141,7 @@ namespace polyhedralGravity {
          * @note Just used for testing purpose
          */
         bool operator==(const TranscendentalExpression &rhs) const {
-            return util::floatNear(ln, rhs.ln) && util::floatNear(an, rhs.an);
+            return util::almostEqualRelative(ln, rhs.ln) && util::almostEqualRelative(an, rhs.an);
         }
 
         /**
@@ -202,10 +202,10 @@ namespace polyhedralGravity {
          * @note Just used for testing purpose
          */
         bool operator==(const HessianPlane &rhs) const {
-            return util::floatNear(a, rhs.a) &&
-                   util::floatNear(b, rhs.b) &&
-                   util::floatNear(c, rhs.c) &&
-                   util::floatNear(d, rhs.d);
+            return util::almostEqualRelative(a, rhs.a) &&
+                   util::almostEqualRelative(b, rhs.b) &&
+                   util::almostEqualRelative(c, rhs.c) &&
+                   util::almostEqualRelative(d, rhs.d);
         }
 
         /**

src/polyhedralGravity/model/GravityModelDetail.cpp

@@ -28,7 +28,7 @@ namespace polyhedralGravity::GravityModel::detail {
         //Calculate N_i * -G_i1 where * is the dot product and then use the inverted sgn
         //We abstain on the double multiplication with -1 in the line above and beyond since two
         //times multiplying with -1 equals no change
-        return sgn(dot(planeUnitNormal, vertex0), util::EPSILON);
+        return sgn(dot(planeUnitNormal, vertex0), util::EPSILON_ZERO_OFFSET);
     }
 
     HessianPlane computeHessianPlane(const Array3 &p, const Array3 &q, const Array3 &r) {
@@ -95,7 +95,7 @@ namespace polyhedralGravity::GravityModel::detail {
                        segmentNormalOrientations.begin(),
                        [&projectionPointOnPlane](const Array3 &segmentUnitNormal, const Array3 &vertex) {
                            using namespace util;
-                           return sgn((dot(segmentUnitNormal, projectionPointOnPlane - vertex)), util::EPSILON) * -1.0;
+                           return sgn((dot(segmentUnitNormal, projectionPointOnPlane - vertex)), util::EPSILON_ZERO_OFFSET) * -1.0;
                        });
         return segmentNormalOrientations;
     }
@@ -190,8 +190,8 @@ namespace polyhedralGravity::GravityModel::detail {
 
                               //4. Option: |s1 - l1| == 0 && |s2 - l2| == 0 Computation point P is located from the beginning on
                               // the direction of a specific segment (P coincides with P' and P'')
-                              if (std::abs(distance.s1 - distance.l1) < EPSILON &&
-                                  std::abs(distance.s2 - distance.l2) < EPSILON) {
+                              if (std::abs(distance.s1 - distance.l1) < EPSILON_ZERO_OFFSET &&
+                                  std::abs(distance.s2 - distance.l2) < EPSILON_ZERO_OFFSET) {
                                   //4. Option - Case 2: P is located on the segment from its right side
                                   // s1 = -|s1|, s2 = -|s2|, l1 = -|l1|, l2 = -|l2|
                                   if (distance.s2 < distance.s1) {
@@ -200,7 +200,7 @@ namespace polyhedralGravity::GravityModel::detail {
                                       distance.l1 *= -1.0;
                                       distance.l2 *= -1.0;
                                       return distance;
-                                  } else if (std::abs(distance.s2 - distance.s1) < EPSILON) {
+                                  } else if (std::abs(distance.s2 - distance.s1) < EPSILON_ZERO_OFFSET) {
                                       //4. Option - Case 1: P is located inside the segment (s2 == s1)
                                       // s1 = -|s1|, s2 = |s2|, l1 = -|l1|, l2 = |l2|
                                       distance.s1 *= -1.0;
@@ -264,9 +264,9 @@ namespace polyhedralGravity::GravityModel::detail {
                               // If sigma_pq == 0 && either of the distances of P' to the two segment endpoints == 0 OR
                               // the 1D and 3D distances are smaller than some EPSILON
                               // then LN_pq can be set to zero
-                              if ((segmentNormalOrientation == 0.0 && (r1Norm < EPSILON || r2Norm < EPSILON)) ||
-                                  (std::abs(distance.s1 + distance.s2) < EPSILON &&
-                                   std::abs(distance.l1 + distance.l2) < EPSILON)) {
+                              if ((segmentNormalOrientation == 0.0 && (r1Norm < EPSILON_ZERO_OFFSET || r2Norm < EPSILON_ZERO_OFFSET)) ||
+                                  (std::abs(distance.s1 + distance.s2) < EPSILON_ZERO_OFFSET &&
+                                   std::abs(distance.l1 + distance.l2) < EPSILON_ZERO_OFFSET)) {
                                   transcendentalExpressionPerSegment.ln = 0.0;
                               } else {
                                   //Implementation of
@@ -277,7 +277,7 @@ namespace polyhedralGravity::GravityModel::detail {
 
                               //Compute AN_pq according to (15)
                               // If h_p == 0 or h_pq == 0 then AN_pq is zero, too (distances are always positive!)
-                              if (planeDistance < EPSILON || segmentDistance < EPSILON) {
+                              if (planeDistance < EPSILON_ZERO_OFFSET || segmentDistance < EPSILON_ZERO_OFFSET) {
                                   transcendentalExpressionPerSegment.an = 0.0;
                               } else {
                                   //Implementation of:
@@ -331,7 +331,7 @@ namespace polyhedralGravity::GravityModel::detail {
             const unsigned int j = thrust::get<2>(tuple);
 
             //segmentNormalOrientation != 0.0
-            if (std::abs(segmentNormalOrientation) > EPSILON) {
+            if (std::abs(segmentNormalOrientation) > EPSILON_ZERO_OFFSET) {
                 return false;
             }
 
@@ -358,14 +358,14 @@ namespace polyhedralGravity::GravityModel::detail {
             j = thrust::get<1>(tuple);
 
             //segmentNormalOrientation != 0.0
-            if (std::abs(segmentNormalOrientation) > EPSILON) {
+            if (std::abs(segmentNormalOrientation) > EPSILON_ZERO_OFFSET) {
                 return false;
             }
 
             r1Norm = projectionPointVertexNorms[(j + 1) % 3];
             r2Norm = projectionPointVertexNorms[j];
             //r1Norm == 0.0 || r2Norm == 0.0
-            return r1Norm < EPSILON || r2Norm < EPSILON;
+            return r1Norm < EPSILON_ZERO_OFFSET || r2Norm < EPSILON_ZERO_OFFSET;
         })) {
             using namespace util;
             //Two segment vectors G_1 and G_2 of this plane

src/polyhedralGravity/model/GravityModelDetail.h

@@ -21,6 +21,7 @@
 #include "polyhedralGravity/util/UtilityConstants.h"
 #include "polyhedralGravity/util/UtilityContainer.h"
 #include "polyhedralGravity/util/UtilityThrust.h"
+#include "polyhedralGravity/util/UtilityFloatArithmetic.h"
 #include "polyhedralGravity/output/Logging.h"
 
 namespace polyhedralGravity::GravityModel::detail {

src/polyhedralGravity/model/Polyhedron.cpp

@@ -202,7 +202,7 @@ namespace polyhedralGravity {
         const Array3 rayVector = normal(segmentVector1, segmentVector2);
 
         // The origin of the array has a slight offset in direction of the normal
-        const Array3 rayOrigin = centroid + (rayVector * EPSILON);
+        const Array3 rayOrigin = centroid + (rayVector * EPSILON_ZERO_OFFSET);
 
         // Count every triangular face which is intersected by the ray
         const auto &[begin, end] = this->transformIterator();
@@ -224,7 +224,7 @@ namespace polyhedralGravity {
         const Array3 edge2 = triangle[2] - triangle[0];
         const Array3 h = cross(rayVector, edge2);
         const double a = dot(edge1, h);
-        if (a > -EPSILON && a < EPSILON) {
+        if (a > -EPSILON_ZERO_OFFSET && a < EPSILON_ZERO_OFFSET) {
             return nullptr;
         }
 
@@ -242,7 +242,7 @@ namespace polyhedralGravity {
         }
 
         const double t = f * dot(edge2, q);
-        if (t > EPSILON) {
+        if (t > EPSILON_ZERO_OFFSET) {
             return std::make_unique<Array3>(rayOrigin + rayVector * t);
         } else {
             return nullptr;

src/polyhedralGravity/model/Polyhedron.h

@@ -22,6 +22,7 @@
 #include "thrust/iterator/transform_iterator.h"
 #include "polyhedralGravity/util/UtilityContainer.h"
 #include "polyhedralGravity/util/UtilityConstants.h"
+#include "polyhedralGravity/util/UtilityFloatArithmetic.h"
 
 namespace polyhedralGravity {
 

src/polyhedralGravity/util/UtilityConstants.h

@@ -2,14 +2,6 @@
 
 namespace polyhedralGravity::util {
 
-    /**
-     * The EPSILON used in the polyhedral gravity model.
-     * @related Used to determine if a floating point number is equal to zero as threshold for rounding errors
-     * @related Used for the sgn() function to determine the sign of a double value. Different compilers
-     * produce different results if no EPSILON is applied for the comparison!
-     */
-    constexpr double EPSILON = 1e-14;
-
     /**
      * PI with enough precision
      */

src/polyhedralGravity/util/UtilityFloatArithmetic.cpp

@@ -3,25 +3,49 @@
 
 namespace polyhedralGravity::util {
 
-    bool floatNear(double lhs, double rhs, int ulpDistance) {
+    template<typename FloatType>
+    bool almostEqualUlps(FloatType lhs, FloatType rhs, int ulpDistance) {
+        static_assert(std::is_same_v<FloatType, float> || std::is_same_v<FloatType, double>,
+                      "Template argument must be FloatType: Either float or double!");
+
+        // In case the floats are equal in their representation, return true
         if (lhs == rhs) {
             return true;
         }
-        if (lhs < 0.0 && rhs > 0.0 || lhs > 0.0 && rhs < 0.0) {
+
+        // In case the signs mismatch, return false
+        if (lhs < static_cast<FloatType>(0.0) && rhs > static_cast<FloatType>(0.0) ||
+            lhs > static_cast<FloatType>(0.0) && rhs < static_cast<FloatType>(0.0)) {
             return false;
-        }
-        return reinterpret_cast<std::int64_t&>(rhs) - reinterpret_cast<std::int64_t&>(lhs) <= ulpDistance;
-    }
+            }
 
-    bool floatNear(float lhs, float rhs, int ulpDistance) {
-        if (lhs == rhs) {
-            return true;
+        if constexpr (std::is_same_v<FloatType, float>) {
+            // In case of float, compute ULP distance by interpreting float as 32-bit integer
+            return reinterpret_cast<std::int32_t&>(rhs) - reinterpret_cast<std::int32_t&>(lhs) <= ulpDistance;
         }
-        if (lhs < 0.0 && rhs > 0.0 || lhs > 0.0 && rhs < 0.0) {
-            return false;
+        else if constexpr (std::is_same_v<FloatType, double>) {
+            // In case of double, compute ULP distance by interpreting double as 64-bit integer
+            return reinterpret_cast<std::int64_t&>(rhs) - reinterpret_cast<std::int64_t&>(lhs) <= ulpDistance;
         }
-        return reinterpret_cast<std::int32_t&>(rhs) - reinterpret_cast<std::int32_t&>(lhs) <= ulpDistance;
+
+        // Due to the static_assert above, this should not happen
+        return false;
+    }
+
+    // Template Instantations for float and double (required since definition is in .cpp file,
+    //      also makes the static assert not strictly necessary)
+    template bool almostEqualUlps<float>(float lhs, float rhs, int ulpDistance);
+    template bool almostEqualUlps<double>(double lhs, double rhs, int ulpDistance);
+
+
+    template<typename FloatType>
+    bool almostEqualRelative(FloatType lhs, FloatType rhs, double epsilon) {
+        const FloatType diff = std::abs(rhs - lhs);
+        const FloatType largerValue = std::max(std::abs(rhs), std::abs(lhs));
+        return diff <= largerValue * epsilon;
     }
 
+    template bool almostEqualRelative<float>(float lhs, float rhs, double epsilon);
+    template bool almostEqualRelative<double>(double lhs, double rhs, double epsilon);
 
 }

src/polyhedralGravity/util/UtilityFloatArithmetic.h

@@ -3,12 +3,32 @@
 #include <algorithm>
 #include <cmath>
 #include <cstdint>
+#include <type_traits>
 
 namespace polyhedralGravity::util {
 
     /**
-     * The maximal allowed ULP distance as computed by {@ref polyhedralGravity::util::floatNear}
-     * utilized for FloatingPoint comparisons in this implementation
+     * The EPSILON used in the polyhedral gravity model to determine a radius around zero/ to use as sligth offset.
+     * @related Used to determine if a floating point number is equal to zero as threshold for rounding errors
+     * @related Used for the sgn() function to determine the sign of a double value. Different compilers
+     * produce different results if no EPSILON is applied for the comparison!
+     */
+    constexpr double EPSILON_ZERO_OFFSET = 1e-14;
+
+    /**
+     * This relative EPSILON is utilized ONLY for testing purposes to compare intermediate values to
+     * Tsoulis' reference implementation Fortran.
+     * It is used in the {@ref polyhedralGravity::util::almostEqualRelative} function.
+     *
+     * @note While in theory no difference at all is observed when compiling this program on Linux using GCC on x86_64,
+     *  the intermediate values change when the program is compiled in different environments.
+     *  Hence, this solves the flakiness of the tests when on different plattforms
+     */
+    constexpr double EPSILON_ALMOST_EQUAL = 1e-10;
+
+    /**
+     * The maximal allowed ULP distance utilized for FloatingPoint comparisons using the
+     * {@ref polyhedralGravity::util::almostEqualUlps} function.
      *
      * @see https://randomascii.wordpress.com/2012/02/25/comparing-floating-point-numbers-2012-edition/
      */
@@ -18,31 +38,36 @@ namespace polyhedralGravity::util {
     /**
      * Function for comparing closeness of two floating point numbers using ULP (Units in the Last Place) method.
      *
+     * @tparam FloatType must be either double or float (ensured by static assertion)
      * @param lhs The left hand side floating point number to compare.
      * @param rhs The right hand side floating point number to compare.
      * @param ulpDistance The maximum acceptable ULP distance between the two floating points
-     *  for which they would be considered near each other. This is optional and by default, it will be MAX_ULP_DISTANCE.
+     *      for which they would be considered near each other. This is optional and by default, it will be MAX_ULP_DISTANCE.
      *
      * @return true if the ULP distance between lhs and rhs is less than or equal to the provided ulpDistance value, otherwise, false.
      *  Returns true if both numbers are exactly the same. Returns false if the signs do not match.
-     *
+     * @example The ULP distance between 3.0 and std::nextafter(3.0, INFINITY) would be 1,
+     *      the ULP distance of 3.0 and std::nextafter(std::nextafter(3.0, INFINITY), INFINITY) would be 2, etc.
      * @see https://randomascii.wordpress.com/2012/02/25/comparing-floating-point-numbers-2012-edition/
      */
-    bool floatNear(double lhs, double rhs, int ulpDistance = MAX_ULP_DISTANCE);
+    template<typename FloatType>
+    bool almostEqualUlps(FloatType lhs, FloatType rhs, int ulpDistance = MAX_ULP_DISTANCE);
 
     /**
-     * Function for comparing closeness of two floating point numbers using ULP (Units in the Last Place) method.
-     *
-     * @param lhs The left hand side floating point number to compare.
-     * @param rhs The right hand side floating point number to compare.
-     * @param ulpDistance The maximum acceptable ULP distance between the two floating points
-     *  for which they would be considered near each other. This is optional and by default, it will be MAX_ULP_DISTANCE.
+     * Function to check if two floating point numbers are relatively equal to each other within a given error range or tolerance.
      *
-     * @return true if the ULP distance between lhs and rhs is less than or equal to the provided ulpDistance value, otherwise, false.
-     *  Returns true if both numbers are exactly the same. Returns false if the signs do not match.
+     * @tparam FloatType must be either double or float (ensured by static assertion)
+     * @param lhs The first floating-point number to be compared.
+     * @param rhs The second floating-point number to be compared.
+     * @param epsilon The tolerance for comparison. Two numbers that are less than epsilon apart are considered equal.
+     *                The default value is `EPSILON`.
      *
+     * @return boolean value - Returns `true` if the absolute difference between `lhs` and `rhs` is less than or equal to
+     *                         the relative error factored by the larger of the magnitude of `lhs` and `rhs`. Otherwise, `false`.
      * @see https://randomascii.wordpress.com/2012/02/25/comparing-floating-point-numbers-2012-edition/
      */
-    bool floatNear(float lhs, float rhs, int ulpDistance = MAX_ULP_DISTANCE);
+    template<typename FloatType>
+    bool almostEqualRelative(FloatType lhs, FloatType rhs, double epsilon = EPSILON_ALMOST_EQUAL);
+
 
 }