512 lines
12 KiB
C++
512 lines
12 KiB
C++
|
// Ceres Solver - A fast non-linear least squares minimizer
|
||
|
// Copyright 2015 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 <algorithm>
|
||
|
#include <cmath>
|
||
|
|
||
|
#include "glog/logging.h"
|
||
|
#include "gtest/gtest.h"
|
||
|
#include "ceres/fpclassify.h"
|
||
|
#include "ceres/stringprintf.h"
|
||
|
#include "ceres/test_util.h"
|
||
|
|
||
|
#define VL VLOG(1)
|
||
|
|
||
|
namespace ceres {
|
||
|
namespace internal {
|
||
|
|
||
|
const double kE = 2.71828182845904523536;
|
||
|
|
||
|
typedef Jet<double, 2> J;
|
||
|
|
||
|
// 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;
|
||
|
}
|
||
|
|
||
|
// On a 32-bit optimized build, the mismatch is about 1.4e-14.
|
||
|
double const kTolerance = 1e-13;
|
||
|
|
||
|
void ExpectJetsClose(const J &x, const J &y) {
|
||
|
ExpectClose(x.a, y.a, kTolerance);
|
||
|
ExpectClose(x.v[0], y.v[0], kTolerance);
|
||
|
ExpectClose(x.v[1], y.v[1], kTolerance);
|
||
|
}
|
||
|
|
||
|
TEST(Jet, Jet) {
|
||
|
// Pick arbitrary values for x and y.
|
||
|
J x = MakeJet(2.3, -2.7, 1e-3);
|
||
|
J y = MakeJet(1.7, 0.5, 1e+2);
|
||
|
|
||
|
VL << "x = " << x;
|
||
|
VL << "y = " << y;
|
||
|
|
||
|
{ // Check that log(exp(x)) == x.
|
||
|
J z = exp(x);
|
||
|
J w = log(z);
|
||
|
VL << "z = " << z;
|
||
|
VL << "w = " << w;
|
||
|
ExpectJetsClose(w, x);
|
||
|
}
|
||
|
|
||
|
{ // Check that (x * y) / x == y.
|
||
|
J z = x * y;
|
||
|
J w = z / x;
|
||
|
VL << "z = " << z;
|
||
|
VL << "w = " << w;
|
||
|
ExpectJetsClose(w, y);
|
||
|
}
|
||
|
|
||
|
{ // Check that sqrt(x * x) == x.
|
||
|
J z = x * x;
|
||
|
J w = sqrt(z);
|
||
|
VL << "z = " << z;
|
||
|
VL << "w = " << w;
|
||
|
ExpectJetsClose(w, x);
|
||
|
}
|
||
|
|
||
|
{ // Check that sqrt(y) * sqrt(y) == y.
|
||
|
J z = sqrt(y);
|
||
|
J w = z * z;
|
||
|
VL << "z = " << z;
|
||
|
VL << "w = " << w;
|
||
|
ExpectJetsClose(w, y);
|
||
|
}
|
||
|
|
||
|
{ // Check that cos(2*x) = cos(x)^2 - sin(x)^2
|
||
|
J z = cos(J(2.0) * x);
|
||
|
J w = cos(x)*cos(x) - sin(x)*sin(x);
|
||
|
VL << "z = " << z;
|
||
|
VL << "w = " << w;
|
||
|
ExpectJetsClose(w, z);
|
||
|
}
|
||
|
|
||
|
{ // Check that sin(2*x) = 2*cos(x)*sin(x)
|
||
|
J z = sin(J(2.0) * x);
|
||
|
J w = J(2.0)*cos(x)*sin(x);
|
||
|
VL << "z = " << z;
|
||
|
VL << "w = " << w;
|
||
|
ExpectJetsClose(w, z);
|
||
|
}
|
||
|
|
||
|
{ // Check that cos(x)*cos(x) + sin(x)*sin(x) = 1
|
||
|
J z = cos(x) * cos(x);
|
||
|
J w = sin(x) * sin(x);
|
||
|
VL << "z = " << z;
|
||
|
VL << "w = " << w;
|
||
|
ExpectJetsClose(z + w, J(1.0));
|
||
|
}
|
||
|
|
||
|
{ // Check that atan2(r*sin(t), r*cos(t)) = t.
|
||
|
J t = MakeJet(0.7, -0.3, +1.5);
|
||
|
J r = MakeJet(2.3, 0.13, -2.4);
|
||
|
VL << "t = " << t;
|
||
|
VL << "r = " << r;
|
||
|
|
||
|
J u = atan2(r * sin(t), r * cos(t));
|
||
|
VL << "u = " << u;
|
||
|
|
||
|
ExpectJetsClose(u, t);
|
||
|
}
|
||
|
|
||
|
{ // Check that tan(x) = sin(x) / cos(x).
|
||
|
J z = tan(x);
|
||
|
J w = sin(x) / cos(x);
|
||
|
VL << "z = " << z;
|
||
|
VL << "w = " << w;
|
||
|
ExpectJetsClose(z, w);
|
||
|
}
|
||
|
|
||
|
{ // Check that tan(atan(x)) = x.
|
||
|
J z = tan(atan(x));
|
||
|
J w = x;
|
||
|
VL << "z = " << z;
|
||
|
VL << "w = " << w;
|
||
|
ExpectJetsClose(z, w);
|
||
|
}
|
||
|
|
||
|
{ // Check that cosh(x)*cosh(x) - sinh(x)*sinh(x) = 1
|
||
|
J z = cosh(x) * cosh(x);
|
||
|
J w = sinh(x) * sinh(x);
|
||
|
VL << "z = " << z;
|
||
|
VL << "w = " << w;
|
||
|
ExpectJetsClose(z - w, J(1.0));
|
||
|
}
|
||
|
|
||
|
{ // Check that tanh(x + y) = (tanh(x) + tanh(y)) / (1 + tanh(x) tanh(y))
|
||
|
J z = tanh(x + y);
|
||
|
J w = (tanh(x) + tanh(y)) / (J(1.0) + tanh(x) * tanh(y));
|
||
|
VL << "z = " << z;
|
||
|
VL << "w = " << w;
|
||
|
ExpectJetsClose(z, w);
|
||
|
}
|
||
|
|
||
|
{ // Check that pow(x, 1) == x.
|
||
|
VL << "x = " << x;
|
||
|
|
||
|
J u = pow(x, 1.);
|
||
|
VL << "u = " << u;
|
||
|
|
||
|
ExpectJetsClose(x, u);
|
||
|
}
|
||
|
|
||
|
{ // Check that pow(x, 1) == x.
|
||
|
J y = MakeJet(1, 0.0, 0.0);
|
||
|
VL << "x = " << x;
|
||
|
VL << "y = " << y;
|
||
|
|
||
|
J u = pow(x, y);
|
||
|
VL << "u = " << u;
|
||
|
|
||
|
ExpectJetsClose(x, u);
|
||
|
}
|
||
|
|
||
|
{ // Check that pow(e, log(x)) == x.
|
||
|
J logx = log(x);
|
||
|
|
||
|
VL << "x = " << x;
|
||
|
VL << "y = " << y;
|
||
|
|
||
|
J u = pow(kE, logx);
|
||
|
VL << "u = " << u;
|
||
|
|
||
|
ExpectJetsClose(x, u);
|
||
|
}
|
||
|
|
||
|
{ // Check that pow(e, log(x)) == x.
|
||
|
J logx = log(x);
|
||
|
J e = MakeJet(kE, 0., 0.);
|
||
|
VL << "x = " << x;
|
||
|
VL << "log(x) = " << logx;
|
||
|
|
||
|
J u = pow(e, logx);
|
||
|
VL << "u = " << u;
|
||
|
|
||
|
ExpectJetsClose(x, u);
|
||
|
}
|
||
|
|
||
|
{ // Check that pow(e, log(x)) == x.
|
||
|
J logx = log(x);
|
||
|
J e = MakeJet(kE, 0., 0.);
|
||
|
VL << "x = " << x;
|
||
|
VL << "logx = " << logx;
|
||
|
|
||
|
J u = pow(e, logx);
|
||
|
VL << "u = " << u;
|
||
|
|
||
|
ExpectJetsClose(x, u);
|
||
|
}
|
||
|
|
||
|
{ // Check that pow(x,y) = exp(y*log(x)).
|
||
|
J logx = log(x);
|
||
|
J e = MakeJet(kE, 0., 0.);
|
||
|
VL << "x = " << x;
|
||
|
VL << "logx = " << logx;
|
||
|
|
||
|
J u = pow(e, y*logx);
|
||
|
J v = pow(x, y);
|
||
|
VL << "u = " << u;
|
||
|
VL << "v = " << v;
|
||
|
|
||
|
ExpectJetsClose(v, u);
|
||
|
}
|
||
|
|
||
|
{ // Check that pow(0, y) == 0 for y > 1, with both arguments Jets.
|
||
|
// This tests special case handling inside pow().
|
||
|
J a = MakeJet(0, 1, 2);
|
||
|
J b = MakeJet(2, 3, 4);
|
||
|
VL << "a = " << a;
|
||
|
VL << "b = " << b;
|
||
|
|
||
|
J c = pow(a, b);
|
||
|
VL << "a^b = " << c;
|
||
|
ExpectJetsClose(c, MakeJet(0, 0, 0));
|
||
|
}
|
||
|
|
||
|
{ // Check that pow(0, y) == 0 for y == 1, with both arguments Jets.
|
||
|
// This tests special case handling inside pow().
|
||
|
J a = MakeJet(0, 1, 2);
|
||
|
J b = MakeJet(1, 3, 4);
|
||
|
VL << "a = " << a;
|
||
|
VL << "b = " << b;
|
||
|
|
||
|
J c = pow(a, b);
|
||
|
VL << "a^b = " << c;
|
||
|
ExpectJetsClose(c, MakeJet(0, 1, 2));
|
||
|
}
|
||
|
|
||
|
{ // Check that 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
|
||
|
VL << "a = " << a;
|
||
|
VL << "b = " << b;
|
||
|
|
||
|
J c = pow(a, b);
|
||
|
VL << "a^b = " << c;
|
||
|
EXPECT_EQ(c.a, 0.0);
|
||
|
EXPECT_FALSE(IsFinite(c.v[0]));
|
||
|
EXPECT_FALSE(IsFinite(c.v[1]));
|
||
|
}
|
||
|
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
|
||
|
VL << "a = " << a;
|
||
|
VL << "b = " << b;
|
||
|
|
||
|
J c = pow(a, b);
|
||
|
VL << "a^b = " << c;
|
||
|
EXPECT_FALSE(IsFinite(c.a));
|
||
|
EXPECT_FALSE(IsFinite(c.v[0]));
|
||
|
EXPECT_FALSE(IsFinite(c.v[1]));
|
||
|
}
|
||
|
|
||
|
{
|
||
|
// The special case of 0^0 = 1 defined by the C standard.
|
||
|
J a = MakeJet(0, 1, 2);
|
||
|
J b = MakeJet(0, 3, 4);
|
||
|
VL << "a = " << a;
|
||
|
VL << "b = " << b;
|
||
|
|
||
|
J c = pow(a, b);
|
||
|
VL << "a^b = " << c;
|
||
|
EXPECT_EQ(c.a, 1.0);
|
||
|
EXPECT_FALSE(IsFinite(c.v[0]));
|
||
|
EXPECT_FALSE(IsFinite(c.v[1]));
|
||
|
}
|
||
|
}
|
||
|
|
||
|
{ // Check that pow(<0, b) is correct for integer b.
|
||
|
// This tests special case handling inside pow().
|
||
|
J a = MakeJet(-1.5, 3, 4);
|
||
|
|
||
|
// b integer:
|
||
|
for (int i = -10; i <= 10; i++) {
|
||
|
J b = MakeJet(i, 0, 5);
|
||
|
VL << "a = " << a;
|
||
|
VL << "b = " << b;
|
||
|
|
||
|
J c = pow(a, b);
|
||
|
VL << "a^b = " << c;
|
||
|
ExpectClose(c.a, pow(-1.5, i), kTolerance);
|
||
|
EXPECT_TRUE(IsFinite(c.v[0]));
|
||
|
EXPECT_FALSE(IsFinite(c.v[1]));
|
||
|
ExpectClose(c.v[0], i * pow(-1.5, i - 1) * 3.0, kTolerance);
|
||
|
}
|
||
|
}
|
||
|
|
||
|
{ // Check that pow(<0, b) is correct for noninteger b.
|
||
|
// This tests special case handling inside pow().
|
||
|
J a = MakeJet(-1.5, 3, 4);
|
||
|
J b = MakeJet(-2.5, 0, 5);
|
||
|
VL << "a = " << a;
|
||
|
VL << "b = " << b;
|
||
|
|
||
|
J c = pow(a, b);
|
||
|
VL << "a^b = " << c;
|
||
|
EXPECT_FALSE(IsFinite(c.a));
|
||
|
EXPECT_FALSE(IsFinite(c.v[0]));
|
||
|
EXPECT_FALSE(IsFinite(c.v[1]));
|
||
|
}
|
||
|
|
||
|
{
|
||
|
// Check that pow(0,y) == 0 for y == 2, with the second argument a
|
||
|
// Jet. This tests special case handling inside pow().
|
||
|
double a = 0;
|
||
|
J b = MakeJet(2, 3, 4);
|
||
|
VL << "a = " << a;
|
||
|
VL << "b = " << b;
|
||
|
|
||
|
J c = pow(a, b);
|
||
|
VL << "a^b = " << c;
|
||
|
ExpectJetsClose(c, MakeJet(0, 0, 0));
|
||
|
}
|
||
|
|
||
|
{
|
||
|
// Check that pow(<0,y) is correct for integer y. This tests special case
|
||
|
// handling inside pow().
|
||
|
double a = -1.5;
|
||
|
for (int i = -10; i <= 10; i++) {
|
||
|
J b = MakeJet(i, 3, 0);
|
||
|
VL << "a = " << a;
|
||
|
VL << "b = " << b;
|
||
|
|
||
|
J c = pow(a, b);
|
||
|
VL << "a^b = " << c;
|
||
|
ExpectClose(c.a, pow(-1.5, i), kTolerance);
|
||
|
EXPECT_FALSE(IsFinite(c.v[0]));
|
||
|
EXPECT_TRUE(IsFinite(c.v[1]));
|
||
|
ExpectClose(c.v[1], 0, kTolerance);
|
||
|
}
|
||
|
}
|
||
|
|
||
|
{
|
||
|
// Check that pow(<0,y) is correct for noninteger y. This tests special
|
||
|
// case handling inside pow().
|
||
|
double a = -1.5;
|
||
|
J b = MakeJet(-3.14, 3, 0);
|
||
|
VL << "a = " << a;
|
||
|
VL << "b = " << b;
|
||
|
|
||
|
J c = pow(a, b);
|
||
|
VL << "a^b = " << c;
|
||
|
EXPECT_FALSE(IsFinite(c.a));
|
||
|
EXPECT_FALSE(IsFinite(c.v[0]));
|
||
|
EXPECT_FALSE(IsFinite(c.v[1]));
|
||
|
}
|
||
|
|
||
|
{ // Check that 1 + x == x + 1.
|
||
|
J a = x + 1.0;
|
||
|
J b = 1.0 + x;
|
||
|
|
||
|
ExpectJetsClose(a, b);
|
||
|
}
|
||
|
|
||
|
{ // Check that 1 - x == -(x - 1).
|
||
|
J a = 1.0 - x;
|
||
|
J b = -(x - 1.0);
|
||
|
|
||
|
ExpectJetsClose(a, b);
|
||
|
}
|
||
|
|
||
|
{ // Check that x/s == x*s.
|
||
|
J a = x / 5.0;
|
||
|
J b = x * 5.0;
|
||
|
|
||
|
ExpectJetsClose(5.0 * a, b / 5.0);
|
||
|
}
|
||
|
|
||
|
{ // Check that x / y == 1 / (y / x).
|
||
|
J a = x / y;
|
||
|
J b = 1.0 / (y / x);
|
||
|
VL << "a = " << a;
|
||
|
VL << "b = " << b;
|
||
|
|
||
|
ExpectJetsClose(a, b);
|
||
|
}
|
||
|
|
||
|
{ // Check that abs(-x * x) == sqrt(x * x).
|
||
|
ExpectJetsClose(abs(-x), sqrt(x * x));
|
||
|
}
|
||
|
|
||
|
{ // Check that cos(acos(x)) == x.
|
||
|
J a = MakeJet(0.1, -2.7, 1e-3);
|
||
|
ExpectJetsClose(cos(acos(a)), a);
|
||
|
ExpectJetsClose(acos(cos(a)), a);
|
||
|
|
||
|
J b = MakeJet(0.6, 0.5, 1e+2);
|
||
|
ExpectJetsClose(cos(acos(b)), b);
|
||
|
ExpectJetsClose(acos(cos(b)), b);
|
||
|
}
|
||
|
|
||
|
{ // Check that sin(asin(x)) == x.
|
||
|
J a = MakeJet(0.1, -2.7, 1e-3);
|
||
|
ExpectJetsClose(sin(asin(a)), a);
|
||
|
ExpectJetsClose(asin(sin(a)), a);
|
||
|
|
||
|
J b = MakeJet(0.4, 0.5, 1e+2);
|
||
|
ExpectJetsClose(sin(asin(b)), b);
|
||
|
ExpectJetsClose(asin(sin(b)), b);
|
||
|
}
|
||
|
}
|
||
|
|
||
|
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;
|
||
|
|
||
|
// Check that M * v == (v^T * M^T)^T
|
||
|
r1 = M * v;
|
||
|
r2 = (v.transpose() * M.transpose()).transpose();
|
||
|
|
||
|
ExpectJetsClose(r1(0), r2(0));
|
||
|
ExpectJetsClose(r1(1), r2(1));
|
||
|
}
|
||
|
|
||
|
TEST(JetTraitsTest, ClassificationMixed) {
|
||
|
Jet<double, 3> a(5.5, 0);
|
||
|
a.v[0] = std::numeric_limits<double>::quiet_NaN();
|
||
|
a.v[1] = std::numeric_limits<double>::infinity();
|
||
|
a.v[2] = -std::numeric_limits<double>::infinity();
|
||
|
EXPECT_FALSE(IsFinite(a));
|
||
|
EXPECT_FALSE(IsNormal(a));
|
||
|
EXPECT_TRUE(IsInfinite(a));
|
||
|
EXPECT_TRUE(IsNaN(a));
|
||
|
}
|
||
|
|
||
|
TEST(JetTraitsTest, ClassificationNaN) {
|
||
|
Jet<double, 3> a(5.5, 0);
|
||
|
a.v[0] = std::numeric_limits<double>::quiet_NaN();
|
||
|
a.v[1] = 0.0;
|
||
|
a.v[2] = 0.0;
|
||
|
EXPECT_FALSE(IsFinite(a));
|
||
|
EXPECT_FALSE(IsNormal(a));
|
||
|
EXPECT_FALSE(IsInfinite(a));
|
||
|
EXPECT_TRUE(IsNaN(a));
|
||
|
}
|
||
|
|
||
|
TEST(JetTraitsTest, ClassificationInf) {
|
||
|
Jet<double, 3> a(5.5, 0);
|
||
|
a.v[0] = std::numeric_limits<double>::infinity();
|
||
|
a.v[1] = 0.0;
|
||
|
a.v[2] = 0.0;
|
||
|
EXPECT_FALSE(IsFinite(a));
|
||
|
EXPECT_FALSE(IsNormal(a));
|
||
|
EXPECT_TRUE(IsInfinite(a));
|
||
|
EXPECT_FALSE(IsNaN(a));
|
||
|
}
|
||
|
|
||
|
TEST(JetTraitsTest, ClassificationFinite) {
|
||
|
Jet<double, 3> a(5.5, 0);
|
||
|
a.v[0] = 100.0;
|
||
|
a.v[1] = 1.0;
|
||
|
a.v[2] = 3.14159;
|
||
|
EXPECT_TRUE(IsFinite(a));
|
||
|
EXPECT_TRUE(IsNormal(a));
|
||
|
EXPECT_FALSE(IsInfinite(a));
|
||
|
EXPECT_FALSE(IsNaN(a));
|
||
|
}
|
||
|
|
||
|
} // namespace internal
|
||
|
} // namespace ceres
|