Extend support for float/double tolerances in testing
authorTeemu Murtola <teemu.murtola@gmail.com>
Sun, 31 Aug 2014 03:54:50 +0000 (06:54 +0300)
committerTeemu Murtola <teemu.murtola@gmail.com>
Sat, 13 Sep 2014 03:20:23 +0000 (05:20 +0200)
Comparisons based only on ULP are awkward when comparing quantities
computed in single and double precision when the accuracy of the
algorithm / implementation does not vary significantly with the
precision.  For example, the total energy when using an Ewald method
will vary with precision, but an acceptable result has nothing to do
with the precision in which it is computed.

Add relativeToleranceAsFloatingPoint to permit expressing the tolerance
as a floating-point number to mean a relative tolerance, but the
implementation is still expressed in ULP.  Make the tolerance class
support different tolerances for single and double precision numbers to
allow this, as well as for other flexibility.

Add tests for the various tolerance construction functions.

Chang variables containing ULP values to be gmx_uint64_t.

Change-Id: I25f7ddf8f4f9a5f1317b1cafd0159d8dd7d62b8c

src/gromacs/fft/tests/fft.cpp
src/testutils/testasserts.cpp
src/testutils/testasserts.h
src/testutils/tests/testasserts_tests.cpp

index 446d50307f9e2bf2da800be03ed2a2550f264f2d..dfe4dbd3311a2ad01ac646e8c88457340e7b6141 100644 (file)
@@ -93,11 +93,8 @@ class BaseFFTTest : public ::testing::Test
             // TODO: These tolerances are just something that has been observed
             // to be sufficient to pass the tests.  It would be nicer to
             // actually argue about why they are sufficient (or what is).
-#ifdef GMX_DOUBLE
-            checker_.setDefaultTolerance(gmx::test::relativeRealTolerance(10.0, 512));
-#else
-            checker_.setDefaultTolerance(gmx::test::relativeRealTolerance(10.0, 64));
-#endif
+            checker_.setDefaultTolerance(
+                    gmx::test::relativeToleranceAsPrecisionDependentUlp(10.0, 64, 512));
         }
         ~BaseFFTTest()
         {
index f14bbe1c5af81eb4bc4494ac5b44ee1e6350faf6..68ea5c720d4fb4ba6bab98cf5b8fdd808ba4e196 100644 (file)
@@ -98,8 +98,8 @@ floatingPointToBiasedInteger(const FloatingPoint<FloatType> &value)
 }
 
 /*! \brief
- * Computes difference in ULPs between two numbers, treating also values of
- * different sign.
+ * Computes the magnitude of the difference in ULPs between two numbers,
+ * treating also values of different sign.
  */
 template <typename FloatType>
 gmx_uint64_t calculateUlpDifference(const FloatingPoint<FloatType> &value1,
@@ -134,6 +134,17 @@ void initDifference(FloatType raw1, FloatType raw2, double *absoluteDifference,
     *ulpDifference      = calculateUlpDifference(value1, value2);
 }
 
+/*! \brief
+ * Converts a relative tolerance into an ULP difference.
+ */
+template <typename FloatType>
+gmx_uint64_t relativeToleranceToUlp(FloatType tolerance)
+{
+    FloatingPoint<FloatType> m(1.0);
+    FloatingPoint<FloatType> t(1.0 + tolerance);
+    return calculateUlpDifference<FloatType>(m, t);
+}
+
 //! \}
 //! \}
 
@@ -164,9 +175,11 @@ bool FloatingPointDifference::isNaN() const
 
 std::string FloatingPointDifference::toString() const
 {
-    return formatString("%g (%" GMX_PRIu64 " %s-prec. ULPs)%s",
+    const double eps = isDouble() ? GMX_DOUBLE_EPS : GMX_FLOAT_EPS;
+    return formatString("%g (%" GMX_PRIu64 " %s-prec. ULPs, rel. %.3g)%s",
                         absoluteDifference_, ulpDifference_,
-                        bDouble_ ? "double" : "single",
+                        isDouble() ? "double" : "single",
+                        ulpDifference_ * eps,
                         bSignDifference_ ? ", signs differ" : "");
 }
 
@@ -187,33 +200,43 @@ bool FloatingPointTolerance::isWithin(
         return false;
     }
 
-    if (difference.asAbsolute() < absoluteTolerance_)
+    const double absoluteTolerance
+        = difference.isDouble() ? doubleAbsoluteTolerance_ : singleAbsoluteTolerance_;
+    if (difference.asAbsolute() < absoluteTolerance)
     {
         return true;
     }
 
-    if (ulpTolerance_ >= 0
-        && difference.asUlps() <= static_cast<gmx_uint64_t>(ulpTolerance_))
+    const gmx_uint64_t ulpTolerance
+        = difference.isDouble() ? doubleUlpTolerance_ : singleUlpTolerance_;
+    if (ulpTolerance < GMX_UINT64_MAX && difference.asUlps() <= ulpTolerance)
     {
         return true;
     }
     return false;
 }
 
-std::string FloatingPointTolerance::toString() const
+std::string FloatingPointTolerance::toString(const FloatingPointDifference &difference) const
 {
-    std::string result;
-    if (absoluteTolerance_ > 0.0)
+    std::string        result;
+    const double       absoluteTolerance
+        = difference.isDouble() ? doubleAbsoluteTolerance_ : singleAbsoluteTolerance_;
+    const gmx_uint64_t ulpTolerance
+        = difference.isDouble() ? doubleUlpTolerance_ : singleUlpTolerance_;
+    const double       eps
+        = difference.isDouble() ? GMX_DOUBLE_EPS : GMX_FLOAT_EPS;
+
+    if (absoluteTolerance > 0.0)
     {
-        result.append(formatString("abs. %g", absoluteTolerance_));
+        result.append(formatString("abs. %g", absoluteTolerance));
     }
-    if (ulpTolerance_ >= 0)
+    if (ulpTolerance < GMX_UINT64_MAX)
     {
         if (!result.empty())
         {
             result.append(", ");
         }
-        result.append(formatString("%d ULPs", ulpTolerance_));
+        result.append(formatString("%" GMX_PRIu64 " ULPs (rel. %.3g)", ulpTolerance, ulpTolerance * eps));
     }
     if (bSignMustMatch_)
     {
@@ -226,5 +249,18 @@ std::string FloatingPointTolerance::toString() const
     return result;
 }
 
+// Doxygen does not recognize this as the same function as in the header...
+//! \cond
+FloatingPointTolerance
+relativeToleranceAsFloatingPoint(double magnitude, double tolerance)
+{
+    const double absoluteTolerance = magnitude * tolerance;
+    return FloatingPointTolerance(absoluteTolerance, absoluteTolerance,
+                                  relativeToleranceToUlp<float>(tolerance),
+                                  relativeToleranceToUlp<double>(tolerance),
+                                  false);
+}
+//! \endcond
+
 } // namespace test
 } // namespace gmx
index b4e0a8ac36f03c7111785f3e092704407cc8e2fc..e38600fb8b0ba02e763d80edd7dec9e55de80882 100644 (file)
@@ -40,6 +40,7 @@
  * assertions for:
  *  - exceptions
  *  - floating-point comparison
+ *  - comparison against NULL
  *
  * \if internal
  * \todo
@@ -221,7 +222,11 @@ class FloatingPointDifference
          * 0.0 and -0.0 are treated as positive and negative, respectively.
          */
         bool signsDiffer() const { return bSignDifference_; }
-
+        /*! \brief
+         * Whether the difference is between single- or double-precision
+         * numbers.
+         */
+        bool isDouble() const { return bDouble_; }
         //! Formats the difference as a string for assertion failure messages.
         std::string toString() const;
 
@@ -242,6 +247,10 @@ class FloatingPointDifference
  * Specifies a floating-point comparison tolerance and checks whether a
  * difference is within the tolerance.
  *
+ * The related functions section lists methods that can be construct methods
+ * using less parameters than the full constructor, and with more obvious
+ * semantics.  These should be preferred over using the constructor directly.
+ *
  * Several types of tolerances are possible:
  *  - _absolute tolerance_: difference between the values must be smaller than
  *    the given tolerance for the check to pass.
@@ -250,7 +259,7 @@ class FloatingPointDifference
  *  - _ULP tolerance_: ULP (units of least precision) difference between the
  *    values must be smaller than the given tolerance for the check to pass.
  *    Setting the ULP tolerance to zero requires exact match.
- *    Setting the ULP tolerance to negative disables the ULP check.
+ *    Setting the ULP tolerance to GMX_UINT64_MAX disables the ULP check.
  *    `0.0` and `-0.0` are treated as equal for the ULP check.
  *  - _sign check_: if set, any values that are of different signs fail the
  *    check (note that this also applies to `0.0` and `-0.0`: a value with a
@@ -264,18 +273,19 @@ class FloatingPointDifference
  * check.  In this case, the sign check must succeed for the check to pass,
  * even if other tolerances are satisfied.
  *
- * Currently, the ULP tolerance is not in any particular precision, but the
- * interpretation depends on the compared numbers: if floats are compared, then
- * the ULP tolerance specifies single-precision ULPs, and if doubles are
- * compared, then the same number is interpreted as double-precision ULPs.
- *
- * The related functions section lists methods that can be construct methods
- * using less parameters than the full constructor, and with more obvious
- * semantics.
+ * The tolerances can be specified separately for single and double precision
+ * comparison.  Different initialization functions have different semantics on
+ * how the provided tolerance values are interpreted; check their
+ * documentation.
  *
  * Methods in this class do not throw, except for toString(), which may throw
  * std::bad_alloc.
  *
+ * \todo
+ * The factory methods that take ULP difference could be better formulated as
+ * methods that take the acceptable number of incorrect bits and/or the number
+ * of accurate bits.
+ *
  * \see FloatingPointDifference
  */
 class FloatingPointTolerance
@@ -284,13 +294,26 @@ class FloatingPointTolerance
         /*! \brief
          * Creates a tolerance with the specified values.
          *
-         * \param[in] absolute       Allowed absolute difference.
-         * \param[in] ulp            Allowed ULP difference.
-         * \param[in] bSignMustMatch Whether sign mismatch fails the comparison.
+         * \param[in]  singleAbsoluteTolerance
+         *     Allowed absolute difference in a single-precision number.
+         * \param[in]  doubleAbsoluteTolerance
+         *     Allowed absolute difference in a double-precision number.
+         * \param[in]  singleUlpTolerance
+         *     Allowed ULP difference in a single-precision number.
+         * \param[in]  doubleUlpTolerance
+         *     Allowed ULP difference in a double-precision number.
+         * \param[in]  bSignMustMatch
+         *     Whether sign mismatch fails the comparison.
          */
-        FloatingPointTolerance(double absolute, int ulp,
-                               bool bSignMustMatch)
-            : absoluteTolerance_(absolute), ulpTolerance_(ulp),
+        FloatingPointTolerance(float        singleAbsoluteTolerance,
+                               double       doubleAbsoluteTolerance,
+                               gmx_uint64_t singleUlpTolerance,
+                               gmx_uint64_t doubleUlpTolerance,
+                               bool         bSignMustMatch)
+            : singleAbsoluteTolerance_(singleAbsoluteTolerance),
+              doubleAbsoluteTolerance_(doubleAbsoluteTolerance),
+              singleUlpTolerance_(singleUlpTolerance),
+              doubleUlpTolerance_(doubleUlpTolerance),
               bSignMustMatch_(bSignMustMatch)
         {
         }
@@ -303,23 +326,79 @@ class FloatingPointTolerance
         bool isWithin(const FloatingPointDifference &difference) const;
 
         //! Formats the tolerance as a string for assertion failure messages.
-        std::string toString() const;
+        std::string toString(const FloatingPointDifference &difference) const;
 
     private:
-        double       absoluteTolerance_;
-        int          ulpTolerance_;
+        float        singleAbsoluteTolerance_;
+        double       doubleAbsoluteTolerance_;
+        gmx_uint64_t singleUlpTolerance_;
+        gmx_uint64_t doubleUlpTolerance_;
         bool         bSignMustMatch_;
 };
 
 /*! \brief
  * Creates a tolerance that only allows a specified ULP difference.
  *
+ * The tolerance uses the given ULP value for both precisions, i.e., double
+ * precision will have much stricter tolerance.
+ *
+ * \related FloatingPointTolerance
+ */
+static inline FloatingPointTolerance
+ulpTolerance(gmx_uint64_t ulpDiff)
+{
+    return FloatingPointTolerance(0.0, 0.0, ulpDiff, ulpDiff, false);
+}
+
+/*! \brief
+ * Creates a tolerance that allows a difference in two compared values that is
+ * relative to the given magnitude.
+ *
+ * \param[in] magnitude  Magnitude of the numbers the computation operates in.
+ * \param[in] tolerance  Relative tolerance permitted (e.g. 1e-4).
+ *
+ * In addition to setting an ULP tolerance equivalent to \p tolerance for both
+ * precisions, this sets the absolute tolerance such that values close to zero
+ * (in general, smaller than \p magnitude) do not fail the check if they
+ * differ by less than \p tolerance evaluated at \p magnitude.  This accounts
+ * for potential loss of precision for small values, and should be used when
+ * accuracy of values much less than \p magnitude do not matter for
+ * correctness.
+ *
+ * The ULP tolerance for different precisions will be different to make them
+ * both match \p tolerance.
+ *
+ * \related FloatingPointTolerance
+ */
+FloatingPointTolerance
+    relativeToleranceAsFloatingPoint(double magnitude, double tolerance);
+
+/*! \brief
+ * Creates a tolerance that allows a precision-dependent relative difference in
+ * a complex computation.
+ *
+ * \param[in] magnitude      Magnitude of the numbers the computation operates in.
+ * \param[in] singleUlpDiff  Expected accuracy of single-precision
+ *     computation (in ULPs).
+ * \param[in] doubleUlpDiff  Expected accuracy of double-precision
+ *     computation (in ULPs).
+ *
+ * This works as relativeToleranceAsUlp(), but allows setting the ULP
+ * difference separately for the different precisions.  This supports
+ * cases where the double-precision calculation can acceptably has a higher ULP
+ * difference, but relaxing the single-precision tolerance would lead to an
+ * unnecessarily loose test.
+ *
  * \related FloatingPointTolerance
  */
 static inline FloatingPointTolerance
-ulpTolerance(gmx_int64_t ulpDiff)
+relativeToleranceAsPrecisionDependentUlp(double       magnitude,
+                                         gmx_uint64_t singleUlpDiff,
+                                         gmx_uint64_t doubleUlpDiff)
 {
-    return FloatingPointTolerance(0.0, ulpDiff, false);
+    return FloatingPointTolerance(magnitude*singleUlpDiff*GMX_FLOAT_EPS,
+                                  magnitude*doubleUlpDiff*GMX_DOUBLE_EPS,
+                                  singleUlpDiff, doubleUlpDiff, false);
 }
 
 /*! \brief
@@ -329,19 +408,19 @@ ulpTolerance(gmx_int64_t ulpDiff)
  * \param[in] magnitude  Magnitude of the numbers the computation operates in.
  * \param[in] ulpDiff    Expected accuracy of the computation (in ULPs).
  *
- * In addition to setting the ULP tolerance, this sets the absolute tolerance
- * such that values close to zero (in general, smaller than \p magnitude) do
- * not fail the check if they differ by less than \p ulpDiff evaluated at
- * \p magniture.  This accounts for potential loss of precision for small
- * values, and should be used when accuracy of values much less than
- * \p magniture do not matter for correctness.
+ * In addition to setting the ULP tolerance as ulpTolerance(), this sets the
+ * absolute tolerance such that values close to zero (in general, smaller than
+ * \p magnitude) do not fail the check if they differ by less than \p ulpDiff
+ * evaluated at \p magnitude.  This accounts for potential loss of precision
+ * for small values, and should be used when accuracy of values much less than
+ * \p magnitude do not matter for correctness.
  *
  * \related FloatingPointTolerance
  */
 static inline FloatingPointTolerance
-relativeRealTolerance(double magnitude, gmx_int64_t ulpDiff)
+relativeToleranceAsUlp(double magnitude, gmx_uint64_t ulpDiff)
 {
-    return FloatingPointTolerance(magnitude*ulpDiff*GMX_REAL_EPS, ulpDiff, false);
+    return relativeToleranceAsPrecisionDependentUlp(magnitude, ulpDiff, ulpDiff);
 }
 
 /*! \brief
@@ -351,7 +430,7 @@ relativeRealTolerance(double magnitude, gmx_int64_t ulpDiff)
  */
 static inline FloatingPointTolerance defaultRealTolerance()
 {
-    return relativeRealTolerance(1.0, 4);
+    return relativeToleranceAsUlp(1.0, 4);
 }
 
 /*! \name Assertions for floating-point comparison
@@ -387,7 +466,7 @@ static inline ::testing::AssertionResult assertEqualWithinTolerance(
            << "  Expected: " << expr1 << std::endl
            << "  Which is: " << value1 << std::endl
            << "Difference: " << diff.toString() << std::endl
-           << " Tolerance: " << tolerance.toString();
+           << " Tolerance: " << tolerance.toString(diff);
 }
 //! \endcond
 
index ed0542adfa575326f001943bb1ea032a7d71d701..01e9d6146da17abb74e72d9d68dc94605c583e89 100644 (file)
@@ -69,6 +69,18 @@ using ::gmx::test::FloatingPointDifference;
 TEST(FloatingPointDifferenceTest, HandlesEqualValues)
 {
     FloatingPointDifference diff(1.2, 1.2);
+    EXPECT_TRUE(diff.isDouble());
+    EXPECT_FALSE(diff.isNaN());
+    EXPECT_EQ(0.0, diff.asAbsolute());
+    EXPECT_EQ(0U,  diff.asUlps());
+    EXPECT_FALSE(diff.signsDiffer());
+}
+
+// TODO: Use typed tests to run all the tests for single and double.
+TEST(FloatingPointDifferenceTest, HandlesFloatValues)
+{
+    FloatingPointDifference diff(1.2f, 1.2f);
+    EXPECT_FALSE(diff.isDouble());
     EXPECT_FALSE(diff.isNaN());
     EXPECT_EQ(0.0, diff.asAbsolute());
     EXPECT_EQ(0U,  diff.asUlps());
@@ -139,4 +151,86 @@ TEST(FloatingPointDifferenceTest, HandlesNaN)
     EXPECT_TRUE(diff.isNaN());
 }
 
+TEST(FloatingPointToleranceTest, UlpTolerance)
+{
+    using gmx::test::ulpTolerance;
+
+    FloatingPointDifference fequal(1.0, 1.0);
+    FloatingPointDifference fulp2(1.0f, addUlps(1.0f, 2));
+    EXPECT_TRUE(ulpTolerance(0).isWithin(fequal));
+    EXPECT_FALSE(ulpTolerance(1).isWithin(fulp2));
+    EXPECT_TRUE(ulpTolerance(2).isWithin(fulp2));
+
+    FloatingPointDifference dequal(1.0, 1.0);
+    FloatingPointDifference dulp2(1.0,  addUlps(1.0, 2));
+    FloatingPointDifference dulp2f(1.0,  static_cast<double>(addUlps(1.0f, 2)));
+    EXPECT_TRUE(ulpTolerance(0).isWithin(dequal));
+    EXPECT_TRUE(ulpTolerance(2).isWithin(dulp2));
+    EXPECT_FALSE(ulpTolerance(2).isWithin(dulp2f));
+}
+
+TEST(FloatingPointToleranceTest, RelativeToleranceAsFloatingPoint)
+{
+    using gmx::test::relativeToleranceAsFloatingPoint;
+
+    FloatingPointDifference fequal(1.0f, 1.0f);
+    FloatingPointDifference fulp2(1.0f, addUlps(1.0f, 2));
+    FloatingPointDifference fdiff(1.0f, 1.011f);
+    FloatingPointDifference fsmall(0.1f, 0.111f);
+    FloatingPointDifference fsmall2(0.1f, 0.121f);
+    EXPECT_TRUE(relativeToleranceAsFloatingPoint(1.0, 1e-2).isWithin(fequal));
+    EXPECT_TRUE(relativeToleranceAsFloatingPoint(1.0, 1e-9).isWithin(fequal));
+    EXPECT_TRUE(relativeToleranceAsFloatingPoint(1.0, 1e-2).isWithin(fulp2));
+    EXPECT_FALSE(relativeToleranceAsFloatingPoint(1.0, 1e-9).isWithin(fulp2));
+    EXPECT_FALSE(relativeToleranceAsFloatingPoint(1.0, 1e-2).isWithin(fdiff));
+    EXPECT_TRUE(relativeToleranceAsFloatingPoint(1.0, 2e-2).isWithin(fdiff));
+    EXPECT_TRUE(relativeToleranceAsFloatingPoint(1.0, 2e-2).isWithin(fsmall));
+    EXPECT_FALSE(relativeToleranceAsFloatingPoint(0.1, 2e-2).isWithin(fsmall));
+    EXPECT_FALSE(relativeToleranceAsFloatingPoint(1.0, 2e-2).isWithin(fsmall2));
+
+    FloatingPointDifference dequal(1.0, 1.0);
+    FloatingPointDifference dulp2f(1.0, static_cast<double>(addUlps(1.0f, 2)));
+    FloatingPointDifference ddiff(1.0, 1.011);
+    FloatingPointDifference dsmall(0.1, 0.111);
+    FloatingPointDifference dsmall2(0.1, 0.121);
+    EXPECT_TRUE(relativeToleranceAsFloatingPoint(1.0, 1e-2).isWithin(dequal));
+    EXPECT_TRUE(relativeToleranceAsFloatingPoint(1.0, 1e-9).isWithin(dequal));
+    EXPECT_TRUE(relativeToleranceAsFloatingPoint(1.0, 1e-2).isWithin(dulp2f));
+    EXPECT_FALSE(relativeToleranceAsFloatingPoint(1.0, 1e-9).isWithin(dulp2f));
+    EXPECT_FALSE(relativeToleranceAsFloatingPoint(1.0, 1e-2).isWithin(ddiff));
+    EXPECT_TRUE(relativeToleranceAsFloatingPoint(1.0, 2e-2).isWithin(ddiff));
+    EXPECT_TRUE(relativeToleranceAsFloatingPoint(1.0, 2e-2).isWithin(dsmall));
+    EXPECT_FALSE(relativeToleranceAsFloatingPoint(0.1, 2e-2).isWithin(dsmall));
+    EXPECT_FALSE(relativeToleranceAsFloatingPoint(1.0, 2e-2).isWithin(dsmall2));
+}
+
+TEST(FloatingPointToleranceTest, RelativeToleranceAsUlp)
+{
+    using gmx::test::relativeToleranceAsUlp;
+
+    FloatingPointDifference fequal(1.0f, 1.0f);
+    FloatingPointDifference fulp4(1.0f, addUlps(1.0f, 4));
+    FloatingPointDifference fsmall(0.1f, addUlps(1.0f, 2) - 0.9f);
+    FloatingPointDifference fsmall2(0.1f, addUlps(1.0f, 6) - 0.9f);
+    EXPECT_TRUE(relativeToleranceAsUlp(1.0, 2).isWithin(fequal));
+    EXPECT_FALSE(relativeToleranceAsUlp(1.0, 2).isWithin(fulp4));
+    EXPECT_TRUE(relativeToleranceAsUlp(1.0, 4).isWithin(fulp4));
+    EXPECT_TRUE(relativeToleranceAsUlp(1.0, 4).isWithin(fsmall));
+    EXPECT_FALSE(relativeToleranceAsUlp(0.1, 4).isWithin(fsmall));
+    EXPECT_FALSE(relativeToleranceAsUlp(1.0, 4).isWithin(fsmall2));
+
+    FloatingPointDifference dequal(1.0, 1.0);
+    FloatingPointDifference dulp4(1.0, addUlps(1.0, 4));
+    FloatingPointDifference dulp4f(1.0,  static_cast<double>(addUlps(1.0f, 4)));
+    FloatingPointDifference dsmall(0.1, addUlps(1.0, 2) - 0.9);
+    FloatingPointDifference dsmall2(0.1, addUlps(1.0, 6) - 0.9);
+    EXPECT_TRUE(relativeToleranceAsUlp(1.0, 2).isWithin(dequal));
+    EXPECT_FALSE(relativeToleranceAsUlp(1.0, 2).isWithin(dulp4));
+    EXPECT_TRUE(relativeToleranceAsUlp(1.0, 4).isWithin(dulp4));
+    EXPECT_FALSE(relativeToleranceAsUlp(1.0, 4).isWithin(dulp4f));
+    EXPECT_TRUE(relativeToleranceAsUlp(1.0, 4).isWithin(dsmall));
+    EXPECT_FALSE(relativeToleranceAsUlp(0.1, 4).isWithin(dsmall));
+    EXPECT_FALSE(relativeToleranceAsUlp(1.0, 4).isWithin(dsmall2));
+}
+
 } // namespace