Resolve two old TODOs in jet.h 1. Use using directive to pull functions from the std:: namespace instead of defining a forwarding function. 2. Remove Eigen 2.0 related definitions. Change-Id: If5f24ef740c17bc4db300fa04f4a2e8f809970d2
diff --git a/include/ceres/jet.h b/include/ceres/jet.h index 2106357..3c53088 100644 --- a/include/ceres/jet.h +++ b/include/ceres/jet.h
@@ -435,35 +435,33 @@ // This is necessary because we want to use the same name (e.g. 'sqrt') for // double-valued and Jet-valued functions, but we are not allowed to put // Jet-valued functions inside namespace std. -// -// TODO(keir): Switch to "using". -inline double abs (double x) { return std::abs(x); } -inline double log (double x) { return std::log(x); } -inline double exp (double x) { return std::exp(x); } -inline double sqrt (double x) { return std::sqrt(x); } -inline double cos (double x) { return std::cos(x); } -inline double acos (double x) { return std::acos(x); } -inline double sin (double x) { return std::sin(x); } -inline double asin (double x) { return std::asin(x); } -inline double tan (double x) { return std::tan(x); } -inline double atan (double x) { return std::atan(x); } -inline double sinh (double x) { return std::sinh(x); } -inline double cosh (double x) { return std::cosh(x); } -inline double tanh (double x) { return std::tanh(x); } -inline double floor (double x) { return std::floor(x); } -inline double ceil (double x) { return std::ceil(x); } -inline double pow (double x, double y) { return std::pow(x, y); } -inline double atan2(double y, double x) { return std::atan2(y, x); } -inline double cbrt (double x) { return std::cbrt(x); } -inline double exp2 (double x) { return std::exp2(x); } -inline double log2 (double x) { return std::log2(x); } -inline double hypot(double x, double y) { return std::hypot(x, y); } -inline double fmax(double x, double y) { return std::fmax(x, y); } -inline double fmin(double x, double y) { return std::fmin(x, y); } -inline double isfinite(double x) { return std::isfinite(x); } -inline double isinf(double x) { return std::isinf(x); } -inline double isnan(double x) { return std::isnan(x); } -inline double isnormal(double x) { return std::isnormal(x); } +using std::abs; +using std::acos; +using std::asin; +using std::atan; +using std::atan2; +using std::cbrt; +using std::ceil; +using std::cos; +using std::cosh; +using std::exp; +using std::exp2; +using std::floor; +using std::fmax; +using std::fmin; +using std::hypot; +using std::isfinite; +using std::isinf; +using std::isnan; +using std::isnormal; +using std::log; +using std::log2; +using std::pow; +using std::sin; +using std::sinh; +using std::sqrt; +using std::tan; +using std::tanh; // Legacy names from pre-C++11 days. inline bool IsFinite (double x) { return std::isfinite(x); } @@ -899,38 +897,6 @@ return Jet<T, N>(tmp1, tmp2 * f.v + tmp3 * g.v); } -// Define the helper functions Eigen needs to embed Jet types. -// -// NOTE(keir): machine_epsilon() and precision() are missing, because they don't -// work with nested template types (e.g. where the scalar is itself templated). -// Among other things, this means that decompositions of Jet's does not work, -// for example -// -// Matrix<Jet<T, N> ... > A, x, b; -// ... -// A.solve(b, &x) -// -// does not work and will fail with a strange compiler error. -// -// TODO(keir): This is an Eigen 2.0 limitation that is lifted in 3.0. When we -// switch to 3.0, also add the rest of the specialization functionality. -template<typename T, int N> inline const Jet<T, N>& ei_conj(const Jet<T, N>& x) { return x; } // NOLINT -template<typename T, int N> inline const Jet<T, N>& ei_real(const Jet<T, N>& x) { return x; } // NOLINT -template<typename T, int N> inline Jet<T, N> ei_imag(const Jet<T, N>& ) { return Jet<T, N>(0.0); } // NOLINT -template<typename T, int N> inline Jet<T, N> ei_abs (const Jet<T, N>& x) { return fabs(x); } // NOLINT -template<typename T, int N> inline Jet<T, N> ei_abs2(const Jet<T, N>& x) { return x * x; } // NOLINT -template<typename T, int N> inline Jet<T, N> ei_sqrt(const Jet<T, N>& x) { return sqrt(x); } // NOLINT -template<typename T, int N> inline Jet<T, N> ei_exp (const Jet<T, N>& x) { return exp(x); } // NOLINT -template<typename T, int N> inline Jet<T, N> ei_log (const Jet<T, N>& x) { return log(x); } // NOLINT -template<typename T, int N> inline Jet<T, N> ei_sin (const Jet<T, N>& x) { return sin(x); } // NOLINT -template<typename T, int N> inline Jet<T, N> ei_cos (const Jet<T, N>& x) { return cos(x); } // NOLINT -template<typename T, int N> inline Jet<T, N> ei_tan (const Jet<T, N>& x) { return tan(x); } // NOLINT -template<typename T, int N> inline Jet<T, N> ei_atan(const Jet<T, N>& x) { return atan(x); } // NOLINT -template<typename T, int N> inline Jet<T, N> ei_sinh(const Jet<T, N>& x) { return sinh(x); } // NOLINT -template<typename T, int N> inline Jet<T, N> ei_cosh(const Jet<T, N>& x) { return cosh(x); } // NOLINT -template<typename T, int N> inline Jet<T, N> ei_tanh(const Jet<T, N>& x) { return tanh(x); } // NOLINT -template<typename T, int N> inline Jet<T, N> ei_pow (const Jet<T, N>& x, Jet<T, N> y) { return pow(x, y); } // NOLINT - // Note: This has to be in the ceres namespace for argument dependent lookup to // function correctly. Otherwise statements like CHECK_LE(x, 2.0) fail with // strange compile errors.
diff --git a/internal/ceres/jet_test.cc b/internal/ceres/jet_test.cc index e15fd67..3c6e23c 100644 --- a/internal/ceres/jet_test.cc +++ b/internal/ceres/jet_test.cc
@@ -30,13 +30,14 @@ #include "ceres/jet.h" +#include <Eigen/Dense> #include <algorithm> #include <cmath> -#include "glog/logging.h" -#include "gtest/gtest.h" #include "ceres/stringprintf.h" #include "ceres/test_util.h" +#include "glog/logging.h" +#include "gtest/gtest.h" #define VL VLOG(1) @@ -760,6 +761,61 @@ EXPECT_FALSE(IsNaN(a)); } +// 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 #if EIGEN_VERSION_AT_LEAST(3, 3, 0) TEST(JetTraitsTest, MatrixScalarUnaryOps) {