2

Some of my unit tests have started failing since adapting some code to enable multi-precision. Header file:

#ifndef SCRATCH_UNITTESTBOOST_INCLUDED
#define SCRATCH_UNITTESTBOOST_INCLUDED

#include <boost/multiprecision/cpp_dec_float.hpp>
// typedef double FLOAT;
typedef boost::multiprecision::cpp_dec_float_50 FLOAT;
const FLOAT ONE(FLOAT(1));

struct Rect
{
    Rect(const FLOAT &width, const FLOAT &height) : Width(width), Height(height){};
    FLOAT getArea() const { return Width * Height; }
    FLOAT Width, Height;
};
#endif

Main test file:

#define BOOST_TEST_DYN_LINK
#define BOOST_TEST_MODULE RectTest
#include <boost/test/unit_test.hpp>
#include "SCRATCH_UnitTestBoost.h"
namespace utf = boost::unit_test;

// Failing
BOOST_AUTO_TEST_CASE(AreaTest1)
{
    Rect R(ONE / 2, ONE / 3);
    FLOAT expected_area = (ONE / 2) * (ONE / 3);

    std::cout << std::setprecision(std::numeric_limits<FLOAT>::digits10) << std::showpoint;
    std::cout << "Expected: " << expected_area << std::endl;
    std::cout << "Actual  : " << R.getArea() << std::endl;

    // BOOST_CHECK_EQUAL(expected_area, R.getArea());
    BOOST_TEST(expected_area == R.getArea());
}

// Tolerance has no effect?
BOOST_AUTO_TEST_CASE(AreaTestTol, *utf::tolerance(1e-40))
{
    Rect R(ONE / 2, ONE / 3);
    FLOAT expected_area = (ONE / 2) * (ONE / 3);
    BOOST_TEST(expected_area == R.getArea());
}

// Passing
BOOST_AUTO_TEST_CASE(AreaTest2)
{
    Rect R(ONE / 7, ONE / 2);
    FLOAT expected_area = (ONE / 7) * (ONE / 2);
    BOOST_CHECK_EQUAL(expected_area, R.getArea());
}

Note that when defining FLOAT as the double type, all the tests pass. What confuses me is that when printing the exact expected and actual values (see AreaTest1) we see the same result. But the error reported from BOOST_TEST is:

    error: in "AreaTest1": check expected_area == R.getArea() has failed 
        [0.16666666666666666666666666666666666666666666666666666666666666666666666666666666 != 
         0.16666666666666666666666666666666666666666666666666666666666666666666666672236366]

Compiling with g++ SCRATCH_UnitTestBoost.cpp -o utb.o -lboost_unit_test_framework.

Questions:

  1. Why is the test failing?
  2. Why does the use of tolerance in AreaTestTol not give outputs as documented here?

Related info:

  1. Tolerances with floating point comparison
  2. Gotchas with multiprecision types
algae
  • 407
  • 4
  • 15

1 Answers1

2

Two Issues:

  • where does the difference come from
  • how to apply the epsilon?

Where The Difference Comes From

Boost Multiprecision uses template expressions to defer evaluation.

Also, you're choosing some rational fractions that cannot be exactly represented base-10 (cpp_dec_float uses decimal, so base-10).

This means that when you do

T x = 1/3;
T y = 1/7;

That will actually approximate both fractions inexactly.

Doing this:

T z = 1/3 * 1/7;

Will actually evaluate the right-handside expression template, so instead of calculating the temporaries like x ans y before, the right hand side has a type of:

expression<detail::multiplies, detail::expression<?>, detail::expression<?>, [2 * ...]>

That's shortened from the actual type:

boost::multiprecision::detail::expression<
    boost::multiprecision::detail::multiplies,
    boost::multiprecision::detail::expression<
        boost::multiprecision::detail::divide_immediates,
        boost::multiprecision::number<boost::multiprecision::backends::cpp_dec_float<50u,
            int, void>, (boost::multiprecision::expression_template_option)1>, int,
            void, void>,
    boost::multiprecision::detail::expression<
        boost::multiprecision::detail::divide_immediates,
        boost::multiprecision::number<boost::multiprecision::backends::cpp_dec_float<50u,
            int, void>, (boost::multiprecision::expression_template_option)1>, int,
            void, void>,
    void, void>

Long story short, this is what you want because it saves you work and keeps better accuracy because the expression is is first normalized to 1/(3*7) so 1/21.

This is where your difference comes from in the first place. Fix it by either:

  1. turning off expression templates

    using T = boost::multiprecision::number<
        boost::multiprecision::cpp_dec_float<50>,
        boost::multiprecision::et_off > >;
    
  2. rewriting the expression to be equivalent of your implementation:

    T expected_area = T(ONE / 7) * T(ONE / 2);
    T expected_area = (ONE / 7).eval() * (ONE / 2).eval();
    

Applying The Tolerance

I find it hard to parse the Boost Unit Test docs on this, but here's empirical data:

BOOST_CHECK_EQUAL(expected_area, R.getArea());
T const eps = std::numeric_limits<T>::epsilon();
BOOST_CHECK_CLOSE(expected_area, R.getArea(), eps);
BOOST_TEST(expected_area == R.getArea(), tt::tolerance(eps));

This fails the first, and passes the last two. Indeed, in addition, the following two also fail:

BOOST_CHECK_EQUAL(expected_area, R.getArea());
BOOST_TEST(expected_area == R.getArea());

So it appears that something has to be done before the utf::tolerance decorator takes effect. Testing with native doubles tells me that only BOOST_TEST applies the tolerance implicitly. So dived into the preprocessed expansion:

    ::boost::unit_test::unit_test_log.set_checkpoint(
        ::boost::unit_test::const_string(
            "/home/sehe/Projects/stackoverflow/test.cpp",
            sizeof("/home/sehe/Projects/stackoverflow/test.cpp") - 1),
        static_cast<std::size_t>(42));
    ::boost::test_tools::tt_detail::report_assertion(
        (::boost::test_tools::assertion::seed()->*a == b).evaluate(),
        (::boost::unit_test::lazy_ostream::instance()
         << ::boost::unit_test::const_string("a == b", sizeof("a == b") - 1)),
        ::boost::unit_test::const_string(
            "/home/sehe/Projects/stackoverflow/test.cpp",
            sizeof("/home/sehe/Projects/stackoverflow/test.cpp") - 1),
        static_cast<std::size_t>(42), ::boost::test_tools::tt_detail::CHECK,
        ::boost::test_tools::tt_detail::CHECK_BUILT_ASSERTION, 0);
} while (::boost::test_tools::tt_detail::dummy_cond());

Digging in a lot more, I ran into:

/*!@brief Indicates if a type can be compared using a tolerance scheme
 *
 * This is a metafunction that should evaluate to @c mpl::true_ if the type
 * @c T can be compared using a tolerance based method, typically for floating point
 * types.
 *
 * This metafunction can be specialized further to declare user types that are
 * floating point (eg. boost.multiprecision).
 */
template <typename T>
struct tolerance_based : tolerance_based_delegate<T, !is_array<T>::value && !is_abstract_class_or_function<T>::value>::type {};

There we have it! But no,

static_assert(boost::math::fpc::tolerance_based<double>::value);
static_assert(boost::math::fpc::tolerance_based<cpp_dec_float_50>::value);

Both already pass. Hmm.

Looking at the decorator I noticed that the tolerance injected into the fixture context is typed.

Experimentally I have reached the conclusion that the tolerance decorator needs to have the same static type argument as the operands in the comparison for it to take effect.

This may actually be very useful (you can have different implicit tolerances for different floating point types), but it is pretty surprising as well.

TL;DR

Here's the full test set fixed and live for your enjoyment:

  • take into account evaluation order and the effect on accuracy
  • use the static type in utf::tolerance(v) to match your operands
  • do not use BOOST_CHECK_EQUAL for tolerance-based comparison
  • I'd suggest to use explicit test_tools::tolerance instead of relying on "ambient" tolerance. After all, we want to be testing our code, not the test framework

Live On Coliru

template <typename T> struct Rect {
    Rect(const T &width, const T &height) : width(width), height(height){};
    T getArea() const { return width * height; }
  private:
    T width, height;
};

#define BOOST_TEST_DYN_LINK
#define BOOST_TEST_MODULE RectTest
#include <boost/multiprecision/cpp_dec_float.hpp>
using DecFloat = boost::multiprecision::cpp_dec_float_50;

#include <boost/test/unit_test.hpp>
namespace utf = boost::unit_test;
namespace tt = boost::test_tools;

namespace {
    template <typename T>
    static inline const T Eps = std::numeric_limits<T>::epsilon();

    template <typename T> struct Fixture {
        T const epsilon = Eps<T>;
        T const ONE     = 1;
        using Rect      = ::Rect<T>;

        void checkArea(int wdenom, int hdenom) const {
            auto w = ONE/wdenom; // could be expression templates
            auto h = ONE/hdenom;

            Rect const R(w, h);
            T expect = w*h;
            BOOST_TEST(expect == R.getArea(), "1/" << wdenom << " x " << "1/" << hdenom);

            // I'd prefer explicit toleranc
            BOOST_TEST(expect == R.getArea(), tt::tolerance(epsilon));
        }
    };

}

BOOST_AUTO_TEST_SUITE(Rectangles)
    BOOST_FIXTURE_TEST_SUITE(Double, Fixture<double>, *utf::tolerance(Eps<double>))
        BOOST_AUTO_TEST_CASE(check2_3)   { checkArea(2, 3); }
        BOOST_AUTO_TEST_CASE(check7_2)   { checkArea(7, 2); }
        BOOST_AUTO_TEST_CASE(check57_31) { checkArea(57, 31); }
    BOOST_AUTO_TEST_SUITE_END()
    BOOST_FIXTURE_TEST_SUITE(MultiPrecision, Fixture<DecFloat>, *utf::tolerance(Eps<DecFloat>))
        BOOST_AUTO_TEST_CASE(check2_3)   { checkArea(2, 3); }
        BOOST_AUTO_TEST_CASE(check7_2)   { checkArea(7, 2); }
        BOOST_AUTO_TEST_CASE(check57_31) { checkArea(57, 31); }
    BOOST_AUTO_TEST_SUITE_END()
BOOST_AUTO_TEST_SUITE_END()

Prints

enter image description here

sehe
  • 374,641
  • 47
  • 450
  • 633
  • What an answer! Not sure I understand the purpose of use cases of expression templates in the first part. However given the terminology, I was able to find [this](https://www.boost.org/doc/libs/1_73_0/libs/math/doc/html/math_toolkit/high_precision/using_test.html) page. – algae Jul 05 '20 at 00:18
  • One other thing.. For an application that switches between multprecision and regular types, should all classes be template classes? Isn't it the same to just have `typedefs` which globally specify the precision type? After all, mixing different objects defined for different types sounds risky. – algae Jul 05 '20 at 00:22
  • 1
    It's the same, but this way I didn't need to recompile to run the different tests. In the interest of test coverage I think that's a good thing to have (because tests that require extra steps can be forgotten). – sehe Jul 05 '20 at 00:36
  • 1
    Also, having templates doesn't mean you need to mix anytthing. Like having `std::vector` as well as `std::vector `doesn't lead to issues. Both can be used separately. I imagine your code-base would have `using Rect = Impl::Rect;` (except I would reserve all-caps for macros) – sehe Jul 05 '20 at 00:36
  • @algae I missed your comment about "expression templates". One of the reasons was in this very post (optimizing evaluation for speed or precision). The page you link tells you how to disable, and that was also in my answer :) If you want an honest understanding of expression templates and how they help: https://www.youtube.com/watch?v=DpxNge0Zhi4 (or start from https://en.wikipedia.org/wiki/Expression_templates) – sehe Jul 05 '20 at 15:50
  • thanks, pity the slides aren't visible. So basically we accumulate results symbolically instead of evaluating operations at each point in code. But making a choice of when to evaluate and how you keep track of lost precision at that point seems non-trivial.. Sounds better to just keep expression templates off and observe the effect of precision program wide, and increase/decrease precision to suit needs – algae Jul 06 '20 at 02:58
  • 1
    You got it. It really just depends on what you're doing. In some cases using expression templates is essential to perf or correctness. I think it only made this particular kind of unit test easier to express exactly, by more predictably sacrificing some precision. – sehe Jul 06 '20 at 03:45