blob: a9207999a2414b45019526aed87930bd07d9236b [file] [log] [blame]
// Ceres Solver - A fast non-linear least squares minimizer
// Copyright 2024 Google Inc. All rights reserved.
// http://ceres-solver.org/
//
// Redistribution and use in source and binary forms, with or without
// modification, are permitted provided that the following conditions are met:
//
// * Redistributions of source code must retain the above copyright notice,
// this list of conditions and the following disclaimer.
// * Redistributions in binary form must reproduce the above copyright notice,
// this list of conditions and the following disclaimer in the documentation
// and/or other materials provided with the distribution.
// * Neither the name of Google Inc. nor the names of its contributors may be
// used to endorse or promote products derived from this software without
// specific prior written permission.
//
// THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS "AS IS"
// AND ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE
// IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE
// ARE DISCLAIMED. IN NO EVENT SHALL THE COPYRIGHT OWNER OR CONTRIBUTORS BE
// LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL, EXEMPLARY, OR
// CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT LIMITED TO, PROCUREMENT OF
// SUBSTITUTE GOODS OR SERVICES; LOSS OF USE, DATA, OR PROFITS; OR BUSINESS
// INTERRUPTION) HOWEVER CAUSED AND ON ANY THEORY OF LIABILITY, WHETHER IN
// CONTRACT, STRICT LIABILITY, OR TORT (INCLUDING NEGLIGENCE OR OTHERWISE)
// ARISING IN ANY WAY OUT OF THE USE OF THIS SOFTWARE, EVEN IF ADVISED OF THE
// POSSIBILITY OF SUCH DAMAGE.
//
// Author: keir@google.com (Keir Mierle)
#include "ceres/jet.h"
#include <Eigen/Dense>
#include <algorithm>
#include <cfenv>
#include <cmath>
#include <stdexcept>
#include "ceres/stringprintf.h"
#include "ceres/test_util.h"
#include "glog/logging.h"
#include "gmock/gmock.h"
#include "gtest/gtest.h"
// The floating-point environment access and modification is only meaningful
// with the following pragma.
#ifdef _MSC_VER
#pragma float_control(precise, on, push)
#pragma fenv_access(on)
#elif !(defined(__ARM_ARCH) && __ARM_ARCH >= 8) && !defined(__MINGW32__)
// NOTE: FENV_ACCESS cannot be set to ON when targeting arm(v8) and MinGW
#pragma STDC FENV_ACCESS ON
#else
#define CERES_NO_FENV_ACCESS
#endif
namespace ceres::internal {
namespace {
constexpr double kE = 2.71828182845904523536;
using J = Jet<double, 2>;
// Don't care about the dual part for scalar part categorization and comparison
// tests
template <typename T>
using J0 = Jet<T, 0>;
using J0d = J0<double>;
// Convenient shorthand for making a jet.
J MakeJet(double a, double v0, double v1) {
J z;
z.a = a;
z.v[0] = v0;
z.v[1] = v1;
return z;
}
constexpr double kTolerance = 1e-13;
// Stores the floating-point environment containing active floating-point
// exceptions, rounding mode, etc., and restores it upon destruction.
//
// Useful for avoiding side-effects.
class Fenv {
public:
Fenv() { std::fegetenv(&e); }
~Fenv() { std::fesetenv(&e); }
Fenv(const Fenv&) = delete;
Fenv& operator=(const Fenv&) = delete;
private:
std::fenv_t e;
};
bool AreAlmostEqual(double x, double y, double max_abs_relative_difference) {
if (std::isnan(x) && std::isnan(y)) {
return true;
}
if (std::isinf(x) && std::isinf(y)) {
return (std::signbit(x) == std::signbit(y));
}
Fenv env; // Do not leak floating-point exceptions to the caller
double absolute_difference = std::abs(x - y);
double relative_difference =
absolute_difference / std::max(std::abs(x), std::abs(y));
if (std::fpclassify(x) == FP_ZERO || std::fpclassify(y) == FP_ZERO) {
// If x or y is exactly zero, then relative difference doesn't have any
// meaning. Take the absolute difference instead.
relative_difference = absolute_difference;
}
return std::islessequal(relative_difference, max_abs_relative_difference);
}
MATCHER_P2(IsAlmostEqualToWithTolerance,
y,
tolerance,
"is almost equal to " + testing::PrintToString(y) +
" with tolerance " + testing::PrintToString(tolerance)) {
const bool result = (AreAlmostEqual(arg.a, y.a, tolerance) &&
AreAlmostEqual(arg.v[0], y.v[0], tolerance) &&
AreAlmostEqual(arg.v[1], y.v[1], tolerance));
if (!result) {
*result_listener << "\nexpected - actual : " << y - arg;
}
return result;
}
MATCHER_P(IsAlmostEqualTo, y, "") {
return ExplainMatchResult(
IsAlmostEqualToWithTolerance(y, kTolerance), arg, result_listener);
}
const double kStep = 1e-8;
const double kNumericalTolerance = 1e-6; // Numeric derivation is quite inexact
// Differentiate using Jet and confirm results with numerical derivation.
template <typename Function>
void NumericalTest(const char* name, const Function& f, const double x) {
const double exact_dx = f(MakeJet(x, 1.0, 0.0)).v[0];
const double estimated_dx =
(f(J(x + kStep)).a - f(J(x - kStep)).a) / (2.0 * kStep);
VLOG(1) << name << "(" << x << "), exact dx: " << exact_dx
<< ", estimated dx: " << estimated_dx;
ExpectClose(exact_dx, estimated_dx, kNumericalTolerance);
}
// Same as NumericalTest, but given a function taking two arguments.
template <typename Function>
void NumericalTest2(const char* name,
const Function& f,
const double x,
const double y) {
const J exact_delta = f(MakeJet(x, 1.0, 0.0), MakeJet(y, 0.0, 1.0));
const double exact_dx = exact_delta.v[0];
const double exact_dy = exact_delta.v[1];
// Sanity check - these should be equivalent:
EXPECT_EQ(exact_dx, f(MakeJet(x, 1.0, 0.0), MakeJet(y, 0.0, 0.0)).v[0]);
EXPECT_EQ(exact_dx, f(MakeJet(x, 0.0, 1.0), MakeJet(y, 0.0, 0.0)).v[1]);
EXPECT_EQ(exact_dy, f(MakeJet(x, 0.0, 0.0), MakeJet(y, 1.0, 0.0)).v[0]);
EXPECT_EQ(exact_dy, f(MakeJet(x, 0.0, 0.0), MakeJet(y, 0.0, 1.0)).v[1]);
const double estimated_dx =
(f(J(x + kStep), J(y)).a - f(J(x - kStep), J(y)).a) / (2.0 * kStep);
const double estimated_dy =
(f(J(x), J(y + kStep)).a - f(J(x), J(y - kStep)).a) / (2.0 * kStep);
VLOG(1) << name << "(" << x << ", " << y << "), exact dx: " << exact_dx
<< ", estimated dx: " << estimated_dx;
ExpectClose(exact_dx, estimated_dx, kNumericalTolerance);
VLOG(1) << name << "(" << x << ", " << y << "), exact dy: " << exact_dy
<< ", estimated dy: " << estimated_dy;
ExpectClose(exact_dy, estimated_dy, kNumericalTolerance);
}
} // namespace
// Pick arbitrary values for x and y.
const J x = MakeJet(2.3, -2.7, 1e-3);
const J y = MakeJet(1.7, 0.5, 1e+2);
const J z = MakeJet(1e-6, 1e-4, 1e-2);
TEST(Jet, Elementary) {
EXPECT_THAT((x * y) / x, IsAlmostEqualTo(y));
EXPECT_THAT(sqrt(x * x), IsAlmostEqualTo(x));
EXPECT_THAT(sqrt(y) * sqrt(y), IsAlmostEqualTo(y));
NumericalTest("sqrt", sqrt<double, 2>, 0.00001);
NumericalTest("sqrt", sqrt<double, 2>, 1.0);
EXPECT_THAT(x + 1.0, IsAlmostEqualTo(1.0 + x));
{
J c = x;
c += 1.0;
EXPECT_THAT(c, IsAlmostEqualTo(1.0 + x));
}
EXPECT_THAT(-(x - 1.0), IsAlmostEqualTo(1.0 - x));
{
J c = x;
c -= 1.0;
EXPECT_THAT(c, IsAlmostEqualTo(x - 1.0));
}
EXPECT_THAT((x * 5.0) / 5.0, IsAlmostEqualTo((x / 5.0) * 5.0));
EXPECT_THAT((x * 5.0) / 5.0, IsAlmostEqualTo(x));
EXPECT_THAT((x / 5.0) * 5.0, IsAlmostEqualTo(x));
{
J c = x;
c /= 5.0;
J d = x;
d *= 5.0;
EXPECT_THAT(c, IsAlmostEqualTo(x / 5.0));
EXPECT_THAT(d, IsAlmostEqualTo(5.0 * x));
}
EXPECT_THAT(1.0 / (y / x), IsAlmostEqualTo(x / y));
}
TEST(Jet, Trigonometric) {
EXPECT_THAT(cos(2.0 * x), IsAlmostEqualTo(cos(x) * cos(x) - sin(x) * sin(x)));
EXPECT_THAT(sin(2.0 * x), IsAlmostEqualTo(2.0 * sin(x) * cos(x)));
EXPECT_THAT(sin(x) * sin(x) + cos(x) * cos(x), IsAlmostEqualTo(J(1.0)));
{
J t = MakeJet(0.7, -0.3, +1.5);
J r = MakeJet(2.3, 0.13, -2.4);
EXPECT_THAT(atan2(r * sin(t), r * cos(t)), IsAlmostEqualTo(t));
}
EXPECT_THAT(sin(x) / cos(x), IsAlmostEqualTo(tan(x)));
EXPECT_THAT(tan(atan(x)), IsAlmostEqualTo(x));
{
J a = MakeJet(0.1, -2.7, 1e-3);
EXPECT_THAT(cos(acos(a)), IsAlmostEqualTo(a));
EXPECT_THAT(acos(cos(a)), IsAlmostEqualTo(a));
J b = MakeJet(0.6, 0.5, 1e+2);
EXPECT_THAT(cos(acos(b)), IsAlmostEqualTo(b));
EXPECT_THAT(acos(cos(b)), IsAlmostEqualTo(b));
}
{
J a = MakeJet(0.1, -2.7, 1e-3);
EXPECT_THAT(sin(asin(a)), IsAlmostEqualTo(a));
EXPECT_THAT(asin(sin(a)), IsAlmostEqualTo(a));
J b = MakeJet(0.4, 0.5, 1e+2);
EXPECT_THAT(sin(asin(b)), IsAlmostEqualTo(b));
EXPECT_THAT(asin(sin(b)), IsAlmostEqualTo(b));
}
}
TEST(Jet, Hyperbolic) {
// cosh(x)*cosh(x) - sinh(x)*sinh(x) = 1
EXPECT_THAT(cosh(x) * cosh(x) - sinh(x) * sinh(x), IsAlmostEqualTo(J(1.0)));
// tanh(x + y) = (tanh(x) + tanh(y)) / (1 + tanh(x) tanh(y))
EXPECT_THAT(
tanh(x + y),
IsAlmostEqualTo((tanh(x) + tanh(y)) / (J(1.0) + tanh(x) * tanh(y))));
}
TEST(Jet, Abs) {
EXPECT_THAT(abs(-x * x), IsAlmostEqualTo(x * x));
EXPECT_THAT(abs(-x), IsAlmostEqualTo(sqrt(x * x)));
{
J a = MakeJet(-std::numeric_limits<double>::quiet_NaN(), 2.0, 4.0);
J b = abs(a);
EXPECT_TRUE(std::signbit(b.v[0]));
EXPECT_TRUE(std::signbit(b.v[1]));
}
}
#if defined(CERES_HAS_POSIX_BESSEL_FUNCTIONS) || \
defined(CERES_HAS_CPP17_BESSEL_FUNCTIONS)
TEST(Jet, Bessel) {
J zero = J(0.0);
J z = MakeJet(0.1, -2.7, 1e-3);
#ifdef CERES_HAS_POSIX_BESSEL_FUNCTIONS
EXPECT_THAT(BesselJ0(zero), IsAlmostEqualTo(J(1.0)));
EXPECT_THAT(BesselJ1(zero), IsAlmostEqualTo(zero));
EXPECT_THAT(BesselJn(2, zero), IsAlmostEqualTo(zero));
EXPECT_THAT(BesselJn(3, zero), IsAlmostEqualTo(zero));
EXPECT_THAT(BesselJ0(z), IsAlmostEqualTo(BesselJn(0, z)));
EXPECT_THAT(BesselJ1(z), IsAlmostEqualTo(BesselJn(1, z)));
// See formula http://dlmf.nist.gov/10.6.E1
EXPECT_THAT(BesselJ0(z) + BesselJn(2, z),
IsAlmostEqualTo((2.0 / z) * BesselJ1(z)));
#endif // CERES_HAS_POSIX_BESSEL_FUNCTIONS
#ifdef CERES_HAS_CPP17_BESSEL_FUNCTIONS
EXPECT_THAT(cyl_bessel_j(0, zero), IsAlmostEqualTo(J(1.0)));
EXPECT_THAT(cyl_bessel_j(1, zero), IsAlmostEqualTo(zero));
EXPECT_THAT(cyl_bessel_j(2, zero), IsAlmostEqualTo(zero));
EXPECT_THAT(cyl_bessel_j(3, zero), IsAlmostEqualTo(zero));
EXPECT_THAT(cyl_bessel_j(0, z), IsAlmostEqualTo(BesselJn(0, z)));
EXPECT_THAT(cyl_bessel_j(1, z), IsAlmostEqualTo(BesselJn(1, z)));
// MSVC Bessel functions and their derivatives produce errors slightly above
// kTolerance. Provide an alternative variant with a relaxed threshold.
constexpr double kRelaxedTolerance = 10 * kTolerance;
// See formula http://dlmf.nist.gov/10.6.E1
EXPECT_THAT(cyl_bessel_j(0, z) + cyl_bessel_j(2, z),
IsAlmostEqualToWithTolerance((2.0 / z) * cyl_bessel_j(1, z),
kRelaxedTolerance));
// MSVC does not throw an exception on invalid first argument
#ifndef _MSC_VER
EXPECT_THROW(cyl_bessel_j(-1, zero), std::domain_error);
#endif // defined(_MSC_VER)
#endif // defined(CERES_HAS_CPP17_BESSEL_FUNCTIONS)
}
#endif // defined(CERES_HAS_POSIX_BESSEL_FUNCTIONS) ||
// defined(CERES_HAS_CPP17_BESSEL_FUNCTIONS)
TEST(Jet, Floor) {
{ // floor of a positive number works.
J a = MakeJet(0.1, -2.7, 1e-3);
J b = floor(a);
J expected = MakeJet(floor(a.a), 0.0, 0.0);
EXPECT_EQ(expected, b);
}
{ // floor of a negative number works.
J a = MakeJet(-1.1, -2.7, 1e-3);
J b = floor(a);
J expected = MakeJet(floor(a.a), 0.0, 0.0);
EXPECT_EQ(expected, b);
}
{ // floor of a positive number works.
J a = MakeJet(10.123, -2.7, 1e-3);
J b = floor(a);
J expected = MakeJet(floor(a.a), 0.0, 0.0);
EXPECT_EQ(expected, b);
}
}
TEST(Jet, Ceil) {
{ // ceil of a positive number works.
J a = MakeJet(0.1, -2.7, 1e-3);
J b = ceil(a);
J expected = MakeJet(ceil(a.a), 0.0, 0.0);
EXPECT_EQ(expected, b);
}
{ // ceil of a negative number works.
J a = MakeJet(-1.1, -2.7, 1e-3);
J b = ceil(a);
J expected = MakeJet(ceil(a.a), 0.0, 0.0);
EXPECT_EQ(expected, b);
}
{ // ceil of a positive number works.
J a = MakeJet(10.123, -2.7, 1e-3);
J b = ceil(a);
J expected = MakeJet(ceil(a.a), 0.0, 0.0);
EXPECT_EQ(expected, b);
}
}
TEST(Jet, Erf) {
{ // erf works.
J a = MakeJet(10.123, -2.7, 1e-3);
J b = erf(a);
J expected = MakeJet(erf(a.a), 0.0, 0.0);
EXPECT_EQ(expected, b);
}
NumericalTest("erf", erf<double, 2>, -1.0);
NumericalTest("erf", erf<double, 2>, 1e-5);
NumericalTest("erf", erf<double, 2>, 0.5);
NumericalTest("erf", erf<double, 2>, 100.0);
}
TEST(Jet, Erfc) {
{ // erfc works.
J a = MakeJet(10.123, -2.7, 1e-3);
J b = erfc(a);
J expected = MakeJet(erfc(a.a), 0.0, 0.0);
EXPECT_EQ(expected, b);
}
NumericalTest("erfc", erfc<double, 2>, -1.0);
NumericalTest("erfc", erfc<double, 2>, 1e-5);
NumericalTest("erfc", erfc<double, 2>, 0.5);
NumericalTest("erfc", erfc<double, 2>, 100.0);
}
TEST(Jet, Cbrt) {
EXPECT_THAT(cbrt(x * x * x), IsAlmostEqualTo(x));
EXPECT_THAT(cbrt(y) * cbrt(y) * cbrt(y), IsAlmostEqualTo(y));
EXPECT_THAT(cbrt(x), IsAlmostEqualTo(pow(x, 1.0 / 3.0)));
NumericalTest("cbrt", cbrt<double, 2>, -1.0);
NumericalTest("cbrt", cbrt<double, 2>, -1e-5);
NumericalTest("cbrt", cbrt<double, 2>, 1e-5);
NumericalTest("cbrt", cbrt<double, 2>, 1.0);
}
TEST(Jet, Log1p) {
EXPECT_THAT(log1p(expm1(x)), IsAlmostEqualTo(x));
EXPECT_THAT(log1p(x), IsAlmostEqualTo(log(J{1} + x)));
{ // log1p(x) does not loose precision for small x
J x = MakeJet(1e-16, 1e-8, 1e-4);
EXPECT_THAT(log1p(x),
IsAlmostEqualTo(MakeJet(9.9999999999999998e-17, 1e-8, 1e-4)));
// log(1 + x) collapses to 0
J v = log(J{1} + x);
EXPECT_TRUE(v.a == 0);
}
}
TEST(Jet, Expm1) {
EXPECT_THAT(expm1(log1p(x)), IsAlmostEqualTo(x));
EXPECT_THAT(expm1(x), IsAlmostEqualTo(exp(x) - 1.0));
{ // expm1(x) does not loose precision for small x
J x = MakeJet(9.9999999999999998e-17, 1e-8, 1e-4);
EXPECT_THAT(expm1(x), IsAlmostEqualTo(MakeJet(1e-16, 1e-8, 1e-4)));
// exp(x) - 1 collapses to 0
J v = exp(x) - J{1};
EXPECT_TRUE(v.a == 0);
}
}
TEST(Jet, Exp2) {
EXPECT_THAT(exp2(x), IsAlmostEqualTo(exp(x * log(2.0))));
NumericalTest("exp2", exp2<double, 2>, -1.0);
NumericalTest("exp2", exp2<double, 2>, -1e-5);
NumericalTest("exp2", exp2<double, 2>, -1e-200);
NumericalTest("exp2", exp2<double, 2>, 0.0);
NumericalTest("exp2", exp2<double, 2>, 1e-200);
NumericalTest("exp2", exp2<double, 2>, 1e-5);
NumericalTest("exp2", exp2<double, 2>, 1.0);
}
TEST(Jet, Log) { EXPECT_THAT(log(exp(x)), IsAlmostEqualTo(x)); }
TEST(Jet, Log10) {
EXPECT_THAT(log10(x), IsAlmostEqualTo(log(x) / log(10)));
NumericalTest("log10", log10<double, 2>, 1e-5);
NumericalTest("log10", log10<double, 2>, 1.0);
NumericalTest("log10", log10<double, 2>, 98.76);
}
TEST(Jet, Log2) {
EXPECT_THAT(log2(x), IsAlmostEqualTo(log(x) / log(2)));
NumericalTest("log2", log2<double, 2>, 1e-5);
NumericalTest("log2", log2<double, 2>, 1.0);
NumericalTest("log2", log2<double, 2>, 100.0);
}
TEST(Jet, Norm) {
EXPECT_THAT(norm(x), IsAlmostEqualTo(x * x));
EXPECT_THAT(norm(-x), IsAlmostEqualTo(x * x));
}
TEST(Jet, Pow) {
EXPECT_THAT(pow(x, 1.0), IsAlmostEqualTo(x));
EXPECT_THAT(pow(x, MakeJet(1.0, 0.0, 0.0)), IsAlmostEqualTo(x));
EXPECT_THAT(pow(kE, log(x)), IsAlmostEqualTo(x));
EXPECT_THAT(pow(MakeJet(kE, 0., 0.), log(x)), IsAlmostEqualTo(x));
EXPECT_THAT(pow(x, y),
IsAlmostEqualTo(pow(MakeJet(kE, 0.0, 0.0), y * log(x))));
// Specially cases
// pow(0, y) == 0 for y > 1, with both arguments Jets.
EXPECT_THAT(pow(MakeJet(0, 1, 2), MakeJet(2, 3, 4)),
IsAlmostEqualTo(MakeJet(0, 0, 0)));
// pow(0, y) == 0 for y == 1, with both arguments Jets.
EXPECT_THAT(pow(MakeJet(0, 1, 2), MakeJet(1, 3, 4)),
IsAlmostEqualTo(MakeJet(0, 1, 2)));
// pow(0, <1) is not finite, with both arguments Jets.
{
for (int i = 1; i < 10; i++) {
J a = MakeJet(0, 1, 2);
J b = MakeJet(i * 0.1, 3, 4); // b = 0.1 ... 0.9
J c = pow(a, b);
EXPECT_EQ(c.a, 0.0) << "\na: " << a << "\nb: " << b << "\na^b: " << c;
EXPECT_FALSE(isfinite(c.v[0]))
<< "\na: " << a << "\nb: " << b << "\na^b: " << c;
EXPECT_FALSE(isfinite(c.v[1]))
<< "\na: " << a << "\nb: " << b << "\na^b: " << c;
}
for (int i = -10; i < 0; i++) {
J a = MakeJet(0, 1, 2);
J b = MakeJet(i * 0.1, 3, 4); // b = -1,-0.9 ... -0.1
J c = pow(a, b);
EXPECT_FALSE(isfinite(c.a))
<< "\na: " << a << "\nb: " << b << "\na^b: " << c;
EXPECT_FALSE(isfinite(c.v[0]))
<< "\na: " << a << "\nb: " << b << "\na^b: " << c;
EXPECT_FALSE(isfinite(c.v[1]))
<< "\na: " << a << "\nb: " << b << "\na^b: " << c;
}
// The special case of 0^0 = 1 defined by the C standard.
{
J a = MakeJet(0, 1, 2);
J b = MakeJet(0, 3, 4);
J c = pow(a, b);
EXPECT_EQ(c.a, 1.0) << "\na: " << a << "\nb: " << b << "\na^b: " << c;
EXPECT_FALSE(isfinite(c.v[0]))
<< "\na: " << a << "\nb: " << b << "\na^b: " << c;
EXPECT_FALSE(isfinite(c.v[1]))
<< "\na: " << a << "\nb: " << b << "\na^b: " << c;
}
}
// pow(<0, b) is correct for integer b.
{
J a = MakeJet(-1.5, 3, 4);
// b integer:
for (int i = -10; i <= 10; i++) {
J b = MakeJet(i, 0, 5);
J c = pow(a, b);
EXPECT_TRUE(AreAlmostEqual(c.a, pow(-1.5, i), kTolerance))
<< "\na: " << a << "\nb: " << b << "\na^b: " << c;
EXPECT_TRUE(isfinite(c.v[0]))
<< "\na: " << a << "\nb: " << b << "\na^b: " << c;
EXPECT_FALSE(isfinite(c.v[1]))
<< "\na: " << a << "\nb: " << b << "\na^b: " << c;
EXPECT_TRUE(
AreAlmostEqual(c.v[0], i * pow(-1.5, i - 1) * 3.0, kTolerance))
<< "\na: " << a << "\nb: " << b << "\na^b: " << c;
}
}
// pow(<0, b) is correct for noninteger b.
{
J a = MakeJet(-1.5, 3, 4);
J b = MakeJet(-2.5, 0, 5);
J c = pow(a, b);
EXPECT_FALSE(isfinite(c.a))
<< "\na: " << a << "\nb: " << b << "\na^b: " << c;
EXPECT_FALSE(isfinite(c.v[0]))
<< "\na: " << a << "\nb: " << b << "\na^b: " << c;
EXPECT_FALSE(isfinite(c.v[1]))
<< "\na: " << a << "\nb: " << b << "\na^b: " << c;
}
// pow(0,y) == 0 for y == 2, with the second argument a Jet.
EXPECT_THAT(pow(0.0, MakeJet(2, 3, 4)), IsAlmostEqualTo(MakeJet(0, 0, 0)));
// pow(<0,y) is correct for integer y.
{
double a = -1.5;
for (int i = -10; i <= 10; i++) {
J b = MakeJet(i, 3, 0);
J c = pow(a, b);
ExpectClose(c.a, pow(-1.5, i), kTolerance);
EXPECT_FALSE(isfinite(c.v[0]))
<< "\na: " << a << "\nb: " << b << "\na^b: " << c;
EXPECT_TRUE(isfinite(c.v[1]))
<< "\na: " << a << "\nb: " << b << "\na^b: " << c;
ExpectClose(c.v[1], 0, kTolerance);
}
}
// pow(<0,y) is correct for noninteger y.
{
double a = -1.5;
J b = MakeJet(-3.14, 3, 0);
J c = pow(a, b);
EXPECT_FALSE(isfinite(c.a))
<< "\na: " << a << "\nb: " << b << "\na^b: " << c;
EXPECT_FALSE(isfinite(c.v[0]))
<< "\na: " << a << "\nb: " << b << "\na^b: " << c;
EXPECT_FALSE(isfinite(c.v[1]))
<< "\na: " << a << "\nb: " << b << "\na^b: " << c;
}
}
TEST(Jet, Hypot2) {
// Resolve the ambiguity between two and three argument hypot overloads
using Hypot2 = J(const J&, const J&);
auto* const hypot2 = static_cast<Hypot2*>(&hypot<double, 2>);
// clang-format off
NumericalTest2("hypot2", hypot2, 0.0, 1e-5);
NumericalTest2("hypot2", hypot2, -1e-5, 0.0);
NumericalTest2("hypot2", hypot2, 1e-5, 1e-5);
NumericalTest2("hypot2", hypot2, 0.0, 1.0);
NumericalTest2("hypot2", hypot2, 1e-3, 1.0);
NumericalTest2("hypot2", hypot2, 1e-3, -1.0);
NumericalTest2("hypot2", hypot2, -1e-3, 1.0);
NumericalTest2("hypot2", hypot2, -1e-3, -1.0);
NumericalTest2("hypot2", hypot2, 1.0, 2.0);
// clang-format on
J zero = MakeJet(0.0, 2.0, 3.14);
EXPECT_THAT(hypot(x, y), IsAlmostEqualTo(sqrt(x * x + y * y)));
EXPECT_THAT(hypot(x, x), IsAlmostEqualTo(sqrt(2.0) * abs(x)));
// The derivative is zero tangentially to the circle:
EXPECT_THAT(hypot(MakeJet(2.0, 1.0, 1.0), MakeJet(2.0, 1.0, -1.0)),
IsAlmostEqualTo(MakeJet(sqrt(8.0), std::sqrt(2.0), 0.0)));
EXPECT_THAT(hypot(zero, x), IsAlmostEqualTo(x));
EXPECT_THAT(hypot(y, zero), IsAlmostEqualTo(y));
// hypot(x, 0, 0) == x, even when x * x underflows:
EXPECT_EQ(
std::numeric_limits<double>::min() * std::numeric_limits<double>::min(),
0.0); // Make sure it underflows
J tiny = MakeJet(std::numeric_limits<double>::min(), 2.0, 3.14);
EXPECT_THAT(hypot(tiny, J{0}), IsAlmostEqualTo(tiny));
// hypot(x, 0, 0) == x, even when x * x overflows:
EXPECT_EQ(
std::numeric_limits<double>::max() * std::numeric_limits<double>::max(),
std::numeric_limits<double>::infinity());
J huge = MakeJet(std::numeric_limits<double>::max(), 2.0, 3.14);
EXPECT_THAT(hypot(huge, J{0}), IsAlmostEqualTo(huge));
}
TEST(Jet, Hypot3) {
J zero = MakeJet(0.0, 2.0, 3.14);
// hypot(x, y, z) == sqrt(x^2 + y^2 + z^2)
EXPECT_THAT(hypot(x, y, z), IsAlmostEqualTo(sqrt(x * x + y * y + z * z)));
// hypot(x, x) == sqrt(3) * abs(x)
EXPECT_THAT(hypot(x, x, x), IsAlmostEqualTo(sqrt(3.0) * abs(x)));
// The derivative is zero tangentially to the circle:
EXPECT_THAT(hypot(MakeJet(2.0, 1.0, 1.0),
MakeJet(2.0, 1.0, -1.0),
MakeJet(2.0, -1.0, 0.0)),
IsAlmostEqualTo(MakeJet(sqrt(12.0), 1.0 / std::sqrt(3.0), 0.0)));
EXPECT_THAT(hypot(x, zero, zero), IsAlmostEqualTo(x));
EXPECT_THAT(hypot(zero, y, zero), IsAlmostEqualTo(y));
EXPECT_THAT(hypot(zero, zero, z), IsAlmostEqualTo(z));
EXPECT_THAT(hypot(x, y, z), IsAlmostEqualTo(hypot(hypot(x, y), z)));
EXPECT_THAT(hypot(x, y, z), IsAlmostEqualTo(hypot(x, hypot(y, z))));
// The following two tests are disabled because the three argument hypot is
// broken in the libc++ shipped with CLANG as of January 2022.
#if !defined(_LIBCPP_VERSION)
// hypot(x, 0, 0) == x, even when x * x underflows:
EXPECT_EQ(
std::numeric_limits<double>::min() * std::numeric_limits<double>::min(),
0.0); // Make sure it underflows
J tiny = MakeJet(std::numeric_limits<double>::min(), 2.0, 3.14);
EXPECT_THAT(hypot(tiny, J{0}, J{0}), IsAlmostEqualTo(tiny));
// hypot(x, 0, 0) == x, even when x * x overflows:
EXPECT_EQ(
std::numeric_limits<double>::max() * std::numeric_limits<double>::max(),
std::numeric_limits<double>::infinity());
J huge = MakeJet(std::numeric_limits<double>::max(), 2.0, 3.14);
EXPECT_THAT(hypot(huge, J{0}, J{0}), IsAlmostEqualTo(huge));
#endif
}
#ifdef CERES_HAS_CPP20
TEST(Jet, Lerp) {
EXPECT_THAT(lerp(x, y, J{0}), IsAlmostEqualTo(x));
EXPECT_THAT(lerp(x, y, J{1}), IsAlmostEqualTo(y));
EXPECT_THAT(lerp(x, x, J{1}), IsAlmostEqualTo(x));
EXPECT_THAT(lerp(y, y, J{0}), IsAlmostEqualTo(y));
EXPECT_THAT(lerp(x, y, J{0.5}), IsAlmostEqualTo((x + y) / J{2.0}));
EXPECT_THAT(lerp(x, y, J{2}), IsAlmostEqualTo(J{2.0} * y - x));
EXPECT_THAT(lerp(x, y, J{-2}), IsAlmostEqualTo(J{3.0} * x - J{2} * y));
}
TEST(Jet, Midpoint) {
EXPECT_THAT(midpoint(x, y), IsAlmostEqualTo((x + y) / J{2}));
EXPECT_THAT(midpoint(x, x), IsAlmostEqualTo(x));
{
// midpoint(x, y) = (x + y) / 2 while avoiding overflow
J x = MakeJet(std::numeric_limits<double>::min(), 1, 2);
J y = MakeJet(std::numeric_limits<double>::max(), 3, 4);
EXPECT_THAT(midpoint(x, y), IsAlmostEqualTo(x + (y - x) / J{2}));
}
{
// midpoint(x, x) = x while avoiding overflow
J x = MakeJet(std::numeric_limits<double>::max(),
std::numeric_limits<double>::max(),
std::numeric_limits<double>::max());
EXPECT_THAT(midpoint(x, x), IsAlmostEqualTo(x));
}
{ // midpoint does not overflow for very large values
constexpr double a = 0.75 * std::numeric_limits<double>::max();
J x = MakeJet(a, a, -a);
J y = MakeJet(a, a, a);
EXPECT_THAT(midpoint(x, y), IsAlmostEqualTo(MakeJet(a, a, 0)));
}
}
#endif // defined(CERES_HAS_CPP20)
TEST(Jet, Fma) {
J v = fma(x, y, z);
J w = x * y + z;
EXPECT_THAT(v, IsAlmostEqualTo(w));
}
TEST(Jet, FmaxJetWithJet) {
Fenv env;
// Clear all exceptions to ensure none are set by the following function
// calls.
std::feclearexcept(FE_ALL_EXCEPT);
EXPECT_THAT(fmax(x, y), IsAlmostEqualTo(x));
EXPECT_THAT(fmax(y, x), IsAlmostEqualTo(x));
// Average the Jets on equality (of scalar parts).
const J scalar_part_only_equal_to_x = J(x.a, 2 * x.v);
const J average = (x + scalar_part_only_equal_to_x) * 0.5;
EXPECT_THAT(fmax(x, scalar_part_only_equal_to_x), IsAlmostEqualTo(average));
EXPECT_THAT(fmax(scalar_part_only_equal_to_x, x), IsAlmostEqualTo(average));
// Follow convention of fmax(): treat NANs as missing values.
const J nan_scalar_part(std::numeric_limits<double>::quiet_NaN(), 2 * x.v);
EXPECT_THAT(fmax(x, nan_scalar_part), IsAlmostEqualTo(x));
EXPECT_THAT(fmax(nan_scalar_part, x), IsAlmostEqualTo(x));
#ifndef CERES_NO_FENV_ACCESS
EXPECT_EQ(std::fetestexcept(FE_ALL_EXCEPT & ~FE_INEXACT), 0);
#endif
}
TEST(Jet, FmaxJetWithScalar) {
Fenv env;
// Clear all exceptions to ensure none are set by the following function
// calls.
std::feclearexcept(FE_ALL_EXCEPT);
EXPECT_THAT(fmax(x, y.a), IsAlmostEqualTo(x));
EXPECT_THAT(fmax(y.a, x), IsAlmostEqualTo(x));
EXPECT_THAT(fmax(y, x.a), IsAlmostEqualTo(J{x.a}));
EXPECT_THAT(fmax(x.a, y), IsAlmostEqualTo(J{x.a}));
// Average the Jet and scalar cast to a Jet on equality (of scalar parts).
const J average = (x + J{x.a}) * 0.5;
EXPECT_THAT(fmax(x, x.a), IsAlmostEqualTo(average));
EXPECT_THAT(fmax(x.a, x), IsAlmostEqualTo(average));
// Follow convention of fmax(): treat NANs as missing values.
EXPECT_THAT(fmax(x, std::numeric_limits<double>::quiet_NaN()),
IsAlmostEqualTo(x));
EXPECT_THAT(fmax(std::numeric_limits<double>::quiet_NaN(), x),
IsAlmostEqualTo(x));
const J nan_scalar_part(std::numeric_limits<double>::quiet_NaN(), 2 * x.v);
EXPECT_THAT(fmax(nan_scalar_part, x.a), IsAlmostEqualTo(J{x.a}));
EXPECT_THAT(fmax(x.a, nan_scalar_part), IsAlmostEqualTo(J{x.a}));
#ifndef CERES_NO_FENV_ACCESS
EXPECT_EQ(std::fetestexcept(FE_ALL_EXCEPT & ~FE_INEXACT), 0);
#endif
}
TEST(Jet, FminJetWithJet) {
Fenv env;
// Clear all exceptions to ensure none are set by the following function
// calls.
std::feclearexcept(FE_ALL_EXCEPT);
EXPECT_THAT(fmin(x, y), IsAlmostEqualTo(y));
EXPECT_THAT(fmin(y, x), IsAlmostEqualTo(y));
// Average the Jets on equality (of scalar parts).
const J scalar_part_only_equal_to_x = J(x.a, 2 * x.v);
const J average = (x + scalar_part_only_equal_to_x) * 0.5;
EXPECT_THAT(fmin(x, scalar_part_only_equal_to_x), IsAlmostEqualTo(average));
EXPECT_THAT(fmin(scalar_part_only_equal_to_x, x), IsAlmostEqualTo(average));
// Follow convention of fmin(): treat NANs as missing values.
const J nan_scalar_part(std::numeric_limits<double>::quiet_NaN(), 2 * x.v);
EXPECT_THAT(fmin(x, nan_scalar_part), IsAlmostEqualTo(x));
EXPECT_THAT(fmin(nan_scalar_part, x), IsAlmostEqualTo(x));
#ifndef CERES_NO_FENV_ACCESS
EXPECT_EQ(std::fetestexcept(FE_ALL_EXCEPT & ~FE_INEXACT), 0);
#endif
}
TEST(Jet, FminJetWithScalar) {
Fenv env;
// Clear all exceptions to ensure none are set by the following function
// calls.
std::feclearexcept(FE_ALL_EXCEPT);
EXPECT_THAT(fmin(x, y.a), IsAlmostEqualTo(J{y.a}));
EXPECT_THAT(fmin(y.a, x), IsAlmostEqualTo(J{y.a}));
EXPECT_THAT(fmin(y, x.a), IsAlmostEqualTo(y));
EXPECT_THAT(fmin(x.a, y), IsAlmostEqualTo(y));
// Average the Jet and scalar cast to a Jet on equality (of scalar parts).
const J average = (x + J{x.a}) * 0.5;
EXPECT_THAT(fmin(x, x.a), IsAlmostEqualTo(average));
EXPECT_THAT(fmin(x.a, x), IsAlmostEqualTo(average));
// Follow convention of fmin(): treat NANs as missing values.
EXPECT_THAT(fmin(x, std::numeric_limits<double>::quiet_NaN()),
IsAlmostEqualTo(x));
EXPECT_THAT(fmin(std::numeric_limits<double>::quiet_NaN(), x),
IsAlmostEqualTo(x));
const J nan_scalar_part(std::numeric_limits<double>::quiet_NaN(), 2 * x.v);
EXPECT_THAT(fmin(nan_scalar_part, x.a), IsAlmostEqualTo(J{x.a}));
EXPECT_THAT(fmin(x.a, nan_scalar_part), IsAlmostEqualTo(J{x.a}));
#ifndef CERES_NO_FENV_ACCESS
EXPECT_EQ(std::fetestexcept(FE_ALL_EXCEPT & ~FE_INEXACT), 0);
#endif
}
TEST(Jet, Fdim) {
Fenv env;
// Clear all exceptions to ensure none are set by the following function
// calls.
std::feclearexcept(FE_ALL_EXCEPT);
const J zero{};
const J diff = x - y;
const J diffx = x - J{y.a};
const J diffy = J{x.a} - y;
EXPECT_THAT(fdim(x, y), IsAlmostEqualTo(diff));
EXPECT_THAT(fdim(y, x), IsAlmostEqualTo(zero));
EXPECT_THAT(fdim(x, y.a), IsAlmostEqualTo(diffx));
EXPECT_THAT(fdim(y.a, x), IsAlmostEqualTo(J{zero.a}));
EXPECT_THAT(fdim(x.a, y), IsAlmostEqualTo(diffy));
EXPECT_THAT(fdim(y, x.a), IsAlmostEqualTo(zero));
EXPECT_TRUE(isnan(fdim(x, std::numeric_limits<J>::quiet_NaN())));
EXPECT_TRUE(isnan(fdim(std::numeric_limits<J>::quiet_NaN(), x)));
EXPECT_TRUE(isnan(fdim(x, std::numeric_limits<double>::quiet_NaN())));
EXPECT_TRUE(isnan(fdim(std::numeric_limits<double>::quiet_NaN(), x)));
#ifndef CERES_NO_FENV_ACCESS
EXPECT_EQ(std::fetestexcept(FE_ALL_EXCEPT & ~FE_INEXACT), 0);
#endif
}
TEST(Jet, CopySign) {
{ // copysign(x, +1)
J z = copysign(x, J{+1});
EXPECT_THAT(z, IsAlmostEqualTo(x));
EXPECT_TRUE(isfinite(z.v[0])) << z;
EXPECT_TRUE(isfinite(z.v[1])) << z;
}
{ // copysign(x, -1)
J z = copysign(x, J{-1});
EXPECT_THAT(z, IsAlmostEqualTo(-x));
EXPECT_TRUE(isfinite(z.v[0])) << z;
EXPECT_TRUE(isfinite(z.v[1])) << z;
}
{ // copysign(-x, +1)
J z = copysign(-x, J{+1});
EXPECT_THAT(z, IsAlmostEqualTo(x));
EXPECT_TRUE(isfinite(z.v[0])) << z;
EXPECT_TRUE(isfinite(z.v[1])) << z;
}
{ // copysign(-x, -1)
J z = copysign(-x, J{-1});
EXPECT_THAT(z, IsAlmostEqualTo(-x));
EXPECT_TRUE(isfinite(z.v[0])) << z;
EXPECT_TRUE(isfinite(z.v[1])) << z;
}
{ // copysign(-0, +1)
J z = copysign(MakeJet(-0, 1, 2), J{+1});
EXPECT_THAT(z, IsAlmostEqualTo(MakeJet(+0, 1, 2)));
EXPECT_FALSE(std::signbit(z.a)) << z;
EXPECT_TRUE(isfinite(z.v[0])) << z;
EXPECT_TRUE(isfinite(z.v[1])) << z;
}
{ // copysign(-0, -1)
J z = copysign(MakeJet(-0, 1, 2), J{-1});
EXPECT_THAT(z, IsAlmostEqualTo(MakeJet(-0, -1, -2)));
EXPECT_TRUE(std::signbit(z.a)) << z;
EXPECT_TRUE(isfinite(z.v[0])) << z;
EXPECT_TRUE(isfinite(z.v[1])) << z;
}
{ // copysign(+0, -1)
J z = copysign(MakeJet(+0, 1, 2), J{-1});
EXPECT_THAT(z, IsAlmostEqualTo(MakeJet(-0, -1, -2)));
EXPECT_TRUE(std::signbit(z.a)) << z;
EXPECT_TRUE(isfinite(z.v[0])) << z;
EXPECT_TRUE(isfinite(z.v[1])) << z;
}
{ // copysign(+0, +1)
J z = copysign(MakeJet(+0, 1, 2), J{+1});
EXPECT_THAT(z, IsAlmostEqualTo(MakeJet(+0, 1, 2)));
EXPECT_FALSE(std::signbit(z.a)) << z;
EXPECT_TRUE(isfinite(z.v[0])) << z;
EXPECT_TRUE(isfinite(z.v[1])) << z;
}
{ // copysign(+0, +0)
J z = copysign(MakeJet(+0, 1, 2), J{+0});
EXPECT_FALSE(std::signbit(z.a)) << z;
EXPECT_TRUE(isnan(z.v[0])) << z;
EXPECT_TRUE(isnan(z.v[1])) << z;
}
{ // copysign(+0, -0)
J z = copysign(MakeJet(+0, 1, 2), J{-0});
EXPECT_FALSE(std::signbit(z.a)) << z;
EXPECT_TRUE(isnan(z.v[0])) << z;
EXPECT_TRUE(isnan(z.v[1])) << z;
}
{ // copysign(-0, +0)
J z = copysign(MakeJet(-0, 1, 2), J{+0});
EXPECT_FALSE(std::signbit(z.a)) << z;
EXPECT_TRUE(isnan(z.v[0])) << z;
EXPECT_TRUE(isnan(z.v[1])) << z;
}
{ // copysign(-0, -0)
J z = copysign(MakeJet(-0, 1, 2), J{-0});
EXPECT_FALSE(std::signbit(z.a)) << z;
EXPECT_TRUE(isnan(z.v[0])) << z;
EXPECT_TRUE(isnan(z.v[1])) << z;
}
{ // copysign(1, -nan)
J z = copysign(MakeJet(1, 2, 3),
-J{std::numeric_limits<double>::quiet_NaN()});
EXPECT_TRUE(std::signbit(z.a)) << z;
EXPECT_TRUE(std::signbit(z.v[0])) << z;
EXPECT_TRUE(std::signbit(z.v[1])) << z;
EXPECT_FALSE(isnan(z.v[0])) << z;
EXPECT_FALSE(isnan(z.v[1])) << z;
}
{ // copysign(1, +nan)
J z = copysign(MakeJet(1, 2, 3),
+J{std::numeric_limits<double>::quiet_NaN()});
EXPECT_FALSE(std::signbit(z.a)) << z;
EXPECT_FALSE(std::signbit(z.v[0])) << z;
EXPECT_FALSE(std::signbit(z.v[1])) << z;
EXPECT_FALSE(isnan(z.v[0])) << z;
EXPECT_FALSE(isnan(z.v[1])) << z;
}
}
TEST(Jet, JetsInEigenMatrices) {
J x = MakeJet(2.3, -2.7, 1e-3);
J y = MakeJet(1.7, 0.5, 1e+2);
J z = MakeJet(5.3, -4.7, 1e-3);
J w = MakeJet(9.7, 1.5, 10.1);
Eigen::Matrix<J, 2, 2> M;
Eigen::Matrix<J, 2, 1> v, r1, r2;
M << x, y, z, w;
v << x, z;
// M * v == (v^T * M^T)^T
r1 = M * v;
r2 = (v.transpose() * M.transpose()).transpose();
EXPECT_THAT(r1(0), IsAlmostEqualTo(r2(0)));
EXPECT_THAT(r1(1), IsAlmostEqualTo(r2(1)));
}
TEST(Jet, ScalarComparison) {
Jet<double, 1> zero{0.0};
zero.v << std::numeric_limits<double>::infinity();
Jet<double, 1> one{1.0};
one.v << std::numeric_limits<double>::quiet_NaN();
Jet<double, 1> two{2.0};
two.v << std::numeric_limits<double>::min() / 2;
Jet<double, 1> three{3.0};
auto inf = std::numeric_limits<Jet<double, 1>>::infinity();
auto nan = std::numeric_limits<Jet<double, 1>>::quiet_NaN();
inf.v << 1.2;
nan.v << 3.4;
std::feclearexcept(FE_ALL_EXCEPT);
EXPECT_FALSE(islessgreater(zero, zero));
EXPECT_FALSE(islessgreater(zero, zero.a));
EXPECT_FALSE(islessgreater(zero.a, zero));
EXPECT_TRUE(isgreaterequal(three, three));
EXPECT_TRUE(isgreaterequal(three, three.a));
EXPECT_TRUE(isgreaterequal(three.a, three));
EXPECT_TRUE(isgreater(three, two));
EXPECT_TRUE(isgreater(three, two.a));
EXPECT_TRUE(isgreater(three.a, two));
EXPECT_TRUE(islessequal(one, one));
EXPECT_TRUE(islessequal(one, one.a));
EXPECT_TRUE(islessequal(one.a, one));
EXPECT_TRUE(isless(one, two));
EXPECT_TRUE(isless(one, two.a));
EXPECT_TRUE(isless(one.a, two));
EXPECT_FALSE(isunordered(inf, one));
EXPECT_FALSE(isunordered(inf, one.a));
EXPECT_FALSE(isunordered(inf.a, one));
EXPECT_TRUE(isunordered(nan, two));
EXPECT_TRUE(isunordered(nan, two.a));
EXPECT_TRUE(isunordered(nan.a, two));
EXPECT_TRUE(isunordered(inf, nan));
EXPECT_TRUE(isunordered(inf, nan.a));
EXPECT_TRUE(isunordered(inf.a, nan.a));
EXPECT_EQ(std::fetestexcept(FE_ALL_EXCEPT & ~FE_INEXACT), 0);
}
TEST(Jet, Nested2XScalarComparison) {
Jet<J0d, 1> zero{J0d{0.0}};
zero.v << std::numeric_limits<J0d>::infinity();
Jet<J0d, 1> one{J0d{1.0}};
one.v << std::numeric_limits<J0d>::quiet_NaN();
Jet<J0d, 1> two{J0d{2.0}};
two.v << std::numeric_limits<J0d>::min() / J0d{2};
Jet<J0d, 1> three{J0d{3.0}};
auto inf = std::numeric_limits<Jet<J0d, 1>>::infinity();
auto nan = std::numeric_limits<Jet<J0d, 1>>::quiet_NaN();
inf.v << J0d{1.2};
nan.v << J0d{3.4};
std::feclearexcept(FE_ALL_EXCEPT);
EXPECT_FALSE(islessgreater(zero, zero));
EXPECT_FALSE(islessgreater(zero, zero.a));
EXPECT_FALSE(islessgreater(zero.a, zero));
EXPECT_FALSE(islessgreater(zero, zero.a.a));
EXPECT_FALSE(islessgreater(zero.a.a, zero));
EXPECT_TRUE(isgreaterequal(three, three));
EXPECT_TRUE(isgreaterequal(three, three.a));
EXPECT_TRUE(isgreaterequal(three.a, three));
EXPECT_TRUE(isgreaterequal(three, three.a.a));
EXPECT_TRUE(isgreaterequal(three.a.a, three));
EXPECT_TRUE(isgreater(three, two));
EXPECT_TRUE(isgreater(three, two.a));
EXPECT_TRUE(isgreater(three.a, two));
EXPECT_TRUE(isgreater(three, two.a.a));
EXPECT_TRUE(isgreater(three.a.a, two));
EXPECT_TRUE(islessequal(one, one));
EXPECT_TRUE(islessequal(one, one.a));
EXPECT_TRUE(islessequal(one.a, one));
EXPECT_TRUE(islessequal(one, one.a.a));
EXPECT_TRUE(islessequal(one.a.a, one));
EXPECT_TRUE(isless(one, two));
EXPECT_TRUE(isless(one, two.a));
EXPECT_TRUE(isless(one.a, two));
EXPECT_TRUE(isless(one, two.a.a));
EXPECT_TRUE(isless(one.a.a, two));
EXPECT_FALSE(isunordered(inf, one));
EXPECT_FALSE(isunordered(inf, one.a));
EXPECT_FALSE(isunordered(inf.a, one));
EXPECT_FALSE(isunordered(inf, one.a.a));
EXPECT_FALSE(isunordered(inf.a.a, one));
EXPECT_TRUE(isunordered(nan, two));
EXPECT_TRUE(isunordered(nan, two.a));
EXPECT_TRUE(isunordered(nan.a, two));
EXPECT_TRUE(isunordered(nan, two.a.a));
EXPECT_TRUE(isunordered(nan.a.a, two));
EXPECT_TRUE(isunordered(inf, nan));
EXPECT_TRUE(isunordered(inf, nan.a));
EXPECT_TRUE(isunordered(inf.a, nan));
EXPECT_TRUE(isunordered(inf, nan.a.a));
EXPECT_TRUE(isunordered(inf.a.a, nan));
EXPECT_EQ(std::fetestexcept(FE_ALL_EXCEPT & ~FE_INEXACT), 0);
}
TEST(JetTraitsTest, ClassificationNaN) {
Jet<double, 1> a(std::numeric_limits<double>::quiet_NaN());
a.v << std::numeric_limits<double>::infinity();
EXPECT_EQ(fpclassify(a), FP_NAN);
EXPECT_FALSE(isfinite(a));
EXPECT_FALSE(isinf(a));
EXPECT_FALSE(isnormal(a));
EXPECT_FALSE(signbit(a));
EXPECT_TRUE(isnan(a));
}
TEST(JetTraitsTest, ClassificationInf) {
Jet<double, 1> a(-std::numeric_limits<double>::infinity());
a.v << std::numeric_limits<double>::quiet_NaN();
EXPECT_EQ(fpclassify(a), FP_INFINITE);
EXPECT_FALSE(isfinite(a));
EXPECT_FALSE(isnan(a));
EXPECT_FALSE(isnormal(a));
EXPECT_TRUE(signbit(a));
EXPECT_TRUE(isinf(a));
}
TEST(JetTraitsTest, ClassificationFinite) {
Jet<double, 1> a(-5.5);
a.v << std::numeric_limits<double>::quiet_NaN();
EXPECT_EQ(fpclassify(a), FP_NORMAL);
EXPECT_FALSE(isinf(a));
EXPECT_FALSE(isnan(a));
EXPECT_TRUE(signbit(a));
EXPECT_TRUE(isfinite(a));
EXPECT_TRUE(isnormal(a));
}
TEST(JetTraitsTest, ClassificationScalar) {
EXPECT_EQ(fpclassify(J0d{+0.0}), FP_ZERO);
EXPECT_EQ(fpclassify(J0d{-0.0}), FP_ZERO);
EXPECT_EQ(fpclassify(J0d{1.234}), FP_NORMAL);
EXPECT_EQ(fpclassify(J0d{std::numeric_limits<double>::min() / 2}),
FP_SUBNORMAL);
EXPECT_EQ(fpclassify(J0d{std::numeric_limits<double>::quiet_NaN()}), FP_NAN);
}
TEST(JetTraitsTest, Nested2XClassificationScalar) {
EXPECT_EQ(fpclassify(J0<J0d>{J0d{+0.0}}), FP_ZERO);
EXPECT_EQ(fpclassify(J0<J0d>{J0d{-0.0}}), FP_ZERO);
EXPECT_EQ(fpclassify(J0<J0d>{J0d{1.234}}), FP_NORMAL);
EXPECT_EQ(fpclassify(J0<J0d>{J0d{std::numeric_limits<double>::min() / 2}}),
FP_SUBNORMAL);
EXPECT_EQ(fpclassify(J0<J0d>{J0d{std::numeric_limits<double>::quiet_NaN()}}),
FP_NAN);
}
// The following test ensures that Jets have all the appropriate Eigen
// related traits so that they can be used as part of matrix
// decompositions.
TEST(Jet, FullRankEigenLLTSolve) {
Eigen::Matrix<J, 3, 3> A;
Eigen::Matrix<J, 3, 1> b, x;
for (int i = 0; i < 3; ++i) {
for (int j = 0; j < 3; ++j) {
A(i, j) = MakeJet(0.0, i, j * j);
}
b(i) = MakeJet(i, i, i);
x(i) = MakeJet(0.0, 0.0, 0.0);
A(i, i) = MakeJet(1.0, i, i * i);
}
x = A.llt().solve(b);
for (int i = 0; i < 3; ++i) {
EXPECT_EQ(x(i).a, b(i).a);
}
}
TEST(Jet, FullRankEigenLDLTSolve) {
Eigen::Matrix<J, 3, 3> A;
Eigen::Matrix<J, 3, 1> b, x;
for (int i = 0; i < 3; ++i) {
for (int j = 0; j < 3; ++j) {
A(i, j) = MakeJet(0.0, i, j * j);
}
b(i) = MakeJet(i, i, i);
x(i) = MakeJet(0.0, 0.0, 0.0);
A(i, i) = MakeJet(1.0, i, i * i);
}
x = A.ldlt().solve(b);
for (int i = 0; i < 3; ++i) {
EXPECT_EQ(x(i).a, b(i).a);
}
}
TEST(Jet, FullRankEigenLUSolve) {
Eigen::Matrix<J, 3, 3> A;
Eigen::Matrix<J, 3, 1> b, x;
for (int i = 0; i < 3; ++i) {
for (int j = 0; j < 3; ++j) {
A(i, j) = MakeJet(0.0, i, j * j);
}
b(i) = MakeJet(i, i, i);
x(i) = MakeJet(0.0, 0.0, 0.0);
A(i, i) = MakeJet(1.0, i, i * i);
}
x = A.lu().solve(b);
for (int i = 0; i < 3; ++i) {
EXPECT_EQ(x(i).a, b(i).a);
}
}
// ScalarBinaryOpTraits is only supported on Eigen versions >= 3.3
TEST(JetTraitsTest, MatrixScalarUnaryOps) {
const J x = MakeJet(2.3, -2.7, 1e-3);
const J y = MakeJet(1.7, 0.5, 1e+2);
Eigen::Matrix<J, 2, 1> a;
a << x, y;
const J sum = a.sum();
const J sum2 = a(0) + a(1);
EXPECT_THAT(sum, IsAlmostEqualTo(sum2));
}
TEST(JetTraitsTest, MatrixScalarBinaryOps) {
const J x = MakeJet(2.3, -2.7, 1e-3);
const J y = MakeJet(1.7, 0.5, 1e+2);
const J z = MakeJet(5.3, -4.7, 1e-3);
const J w = MakeJet(9.7, 1.5, 10.1);
Eigen::Matrix<J, 2, 2> M;
Eigen::Vector2d v;
M << x, y, z, w;
v << 0.6, -2.1;
// M * v == M * v.cast<J>().
const Eigen::Matrix<J, 2, 1> r1 = M * v;
const Eigen::Matrix<J, 2, 1> r2 = M * v.cast<J>();
EXPECT_THAT(r1(0), IsAlmostEqualTo(r2(0)));
EXPECT_THAT(r1(1), IsAlmostEqualTo(r2(1)));
// M * a == M * T(a).
const double a = 3.1;
const Eigen::Matrix<J, 2, 2> r3 = M * a;
const Eigen::Matrix<J, 2, 2> r4 = M * J(a);
EXPECT_THAT(r3(0, 0), IsAlmostEqualTo(r4(0, 0)));
EXPECT_THAT(r3(0, 1), IsAlmostEqualTo(r4(0, 1)));
EXPECT_THAT(r3(1, 0), IsAlmostEqualTo(r4(1, 0)));
EXPECT_THAT(r3(1, 1), IsAlmostEqualTo(r4(1, 1)));
}
TEST(JetTraitsTest, ArrayScalarUnaryOps) {
const J x = MakeJet(2.3, -2.7, 1e-3);
const J y = MakeJet(1.7, 0.5, 1e+2);
Eigen::Array<J, 2, 1> a;
a << x, y;
const J sum = a.sum();
const J sum2 = a(0) + a(1);
EXPECT_THAT(sum, sum2);
}
TEST(JetTraitsTest, ArrayScalarBinaryOps) {
const J x = MakeJet(2.3, -2.7, 1e-3);
const J y = MakeJet(1.7, 0.5, 1e+2);
Eigen::Array<J, 2, 1> a;
Eigen::Array2d b;
a << x, y;
b << 0.6, -2.1;
// a * b == a * b.cast<T>()
const Eigen::Array<J, 2, 1> r1 = a * b;
const Eigen::Array<J, 2, 1> r2 = a * b.cast<J>();
EXPECT_THAT(r1(0), r2(0));
EXPECT_THAT(r1(1), r2(1));
// a * c == a * T(c).
const double c = 3.1;
const Eigen::Array<J, 2, 1> r3 = a * c;
const Eigen::Array<J, 2, 1> r4 = a * J(c);
EXPECT_THAT(r3(0), r3(0));
EXPECT_THAT(r4(1), r4(1));
}
TEST(Jet, Nested3X) {
using JJ = Jet<J, 2>;
using JJJ = Jet<JJ, 2>;
JJJ x;
x.a = JJ(J(1, 0), 0);
x.v[0] = JJ(J(1));
JJJ y = x * x * x;
ExpectClose(y.a.a.a, 1, kTolerance);
ExpectClose(y.v[0].a.a, 3., kTolerance);
ExpectClose(y.v[0].v[0].a, 6., kTolerance);
ExpectClose(y.v[0].v[0].v[0], 6., kTolerance);
JJJ e = exp(x);
ExpectClose(e.a.a.a, kE, kTolerance);
ExpectClose(e.v[0].a.a, kE, kTolerance);
ExpectClose(e.v[0].v[0].a, kE, kTolerance);
ExpectClose(e.v[0].v[0].v[0], kE, kTolerance);
}
#if GTEST_HAS_TYPED_TEST
using Types = testing::Types<std::int16_t,
std::uint16_t,
std::int32_t,
std::uint32_t,
std::int64_t,
std::uint64_t,
float,
double,
long double>;
template <typename T>
class JetTest : public testing::Test {};
TYPED_TEST_SUITE(JetTest, Types);
TYPED_TEST(JetTest, Comparison) {
using Scalar = TypeParam;
EXPECT_EQ(J0<Scalar>{0}, J0<Scalar>{0});
EXPECT_GE(J0<Scalar>{3}, J0<Scalar>{3});
EXPECT_GT(J0<Scalar>{3}, J0<Scalar>{2});
EXPECT_LE(J0<Scalar>{1}, J0<Scalar>{1});
EXPECT_LT(J0<Scalar>{1}, J0<Scalar>{2});
EXPECT_NE(J0<Scalar>{1}, J0<Scalar>{2});
}
TYPED_TEST(JetTest, ScalarComparison) {
using Scalar = TypeParam;
EXPECT_EQ(J0d{0.0}, Scalar{0});
EXPECT_GE(J0d{3.0}, Scalar{3});
EXPECT_GT(J0d{3.0}, Scalar{2});
EXPECT_LE(J0d{1.0}, Scalar{1});
EXPECT_LT(J0d{1.0}, Scalar{2});
EXPECT_NE(J0d{1.0}, Scalar{2});
EXPECT_EQ(Scalar{0}, J0d{0.0});
EXPECT_GE(Scalar{1}, J0d{1.0});
EXPECT_GT(Scalar{2}, J0d{1.0});
EXPECT_LE(Scalar{3}, J0d{3.0});
EXPECT_LT(Scalar{2}, J0d{3.0});
EXPECT_NE(Scalar{2}, J0d{1.0});
}
TYPED_TEST(JetTest, Nested2XComparison) {
using Scalar = TypeParam;
EXPECT_EQ(J0<J0d>{J0d{0.0}}, Scalar{0});
EXPECT_GE(J0<J0d>{J0d{3.0}}, Scalar{3});
EXPECT_GT(J0<J0d>{J0d{3.0}}, Scalar{2});
EXPECT_LE(J0<J0d>{J0d{1.0}}, Scalar{1});
EXPECT_LT(J0<J0d>{J0d{1.0}}, Scalar{2});
EXPECT_NE(J0<J0d>{J0d{1.0}}, Scalar{2});
EXPECT_EQ(Scalar{0}, J0<J0d>{J0d{0.0}});
EXPECT_GE(Scalar{1}, J0<J0d>{J0d{1.0}});
EXPECT_GT(Scalar{2}, J0<J0d>{J0d{1.0}});
EXPECT_LE(Scalar{3}, J0<J0d>{J0d{3.0}});
EXPECT_LT(Scalar{2}, J0<J0d>{J0d{3.0}});
EXPECT_NE(Scalar{2}, J0<J0d>{J0d{1.0}});
}
TYPED_TEST(JetTest, Nested3XComparison) {
using Scalar = TypeParam;
EXPECT_EQ(J0<J0<J0d>>{J0<J0d>{J0d{0.0}}}, Scalar{0});
EXPECT_GE(J0<J0<J0d>>{J0<J0d>{J0d{3.0}}}, Scalar{3});
EXPECT_GT(J0<J0<J0d>>{J0<J0d>{J0d{3.0}}}, Scalar{2});
EXPECT_LE(J0<J0<J0d>>{J0<J0d>{J0d{1.0}}}, Scalar{1});
EXPECT_LT(J0<J0<J0d>>{J0<J0d>{J0d{1.0}}}, Scalar{2});
EXPECT_NE(J0<J0<J0d>>{J0<J0d>{J0d{1.0}}}, Scalar{2});
EXPECT_EQ(Scalar{0}, J0<J0<J0d>>{J0<J0d>{J0d{0.0}}});
EXPECT_GE(Scalar{1}, J0<J0<J0d>>{J0<J0d>{J0d{1.0}}});
EXPECT_GT(Scalar{2}, J0<J0<J0d>>{J0<J0d>{J0d{1.0}}});
EXPECT_LE(Scalar{3}, J0<J0<J0d>>{J0<J0d>{J0d{3.0}}});
EXPECT_LT(Scalar{2}, J0<J0<J0d>>{J0<J0d>{J0d{3.0}}});
EXPECT_NE(Scalar{2}, J0<J0<J0d>>{J0<J0d>{J0d{1.0}}});
}
#endif // GTEST_HAS_TYPED_TEST
} // namespace ceres::internal
#ifdef _MSC_VER
#pragma float_control(pop)
#endif