1// Ceres Solver - A fast non-linear least squares minimizer
2// Copyright 2010, 2011, 2012 Google Inc. All rights reserved.
3// http://code.google.com/p/ceres-solver/
4//
5// Redistribution and use in source and binary forms, with or without
6// modification, are permitted provided that the following conditions are met:
7//
8// * Redistributions of source code must retain the above copyright notice,
9//   this list of conditions and the following disclaimer.
10// * Redistributions in binary form must reproduce the above copyright notice,
11//   this list of conditions and the following disclaimer in the documentation
12//   and/or other materials provided with the distribution.
13// * Neither the name of Google Inc. nor the names of its contributors may be
14//   used to endorse or promote products derived from this software without
15//   specific prior written permission.
16//
17// THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS "AS IS"
18// AND ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE
19// IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE
20// ARE DISCLAIMED. IN NO EVENT SHALL THE COPYRIGHT OWNER OR CONTRIBUTORS BE
21// LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL, EXEMPLARY, OR
22// CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT LIMITED TO, PROCUREMENT OF
23// SUBSTITUTE GOODS OR SERVICES; LOSS OF USE, DATA, OR PROFITS; OR BUSINESS
24// INTERRUPTION) HOWEVER CAUSED AND ON ANY THEORY OF LIABILITY, WHETHER IN
25// CONTRACT, STRICT LIABILITY, OR TORT (INCLUDING NEGLIGENCE OR OTHERWISE)
26// ARISING IN ANY WAY OUT OF THE USE OF THIS SOFTWARE, EVEN IF ADVISED OF THE
27// POSSIBILITY OF SUCH DAMAGE.
28//
29// Author: keir@google.com (Keir Mierle)
30
31#include "ceres/jet.h"
32
33#include <algorithm>
34#include <cmath>
35
36#include "glog/logging.h"
37#include "gtest/gtest.h"
38#include "ceres/fpclassify.h"
39#include "ceres/stringprintf.h"
40#include "ceres/test_util.h"
41
42#define VL VLOG(1)
43
44namespace ceres {
45namespace internal {
46
47const double kE = 2.71828182845904523536;
48
49typedef Jet<double, 2> J;
50
51// Convenient shorthand for making a jet.
52J MakeJet(double a, double v0, double v1) {
53  J z;
54  z.a = a;
55  z.v[0] = v0;
56  z.v[1] = v1;
57  return z;
58}
59
60// On a 32-bit optimized build, the mismatch is about 1.4e-14.
61double const kTolerance = 1e-13;
62
63void ExpectJetsClose(const J &x, const J &y) {
64  ExpectClose(x.a, y.a, kTolerance);
65  ExpectClose(x.v[0], y.v[0], kTolerance);
66  ExpectClose(x.v[1], y.v[1], kTolerance);
67}
68
69TEST(Jet, Jet) {
70  // Pick arbitrary values for x and y.
71  J x = MakeJet(2.3, -2.7, 1e-3);
72  J y = MakeJet(1.7,  0.5, 1e+2);
73
74  VL << "x = " << x;
75  VL << "y = " << y;
76
77  { // Check that log(exp(x)) == x.
78    J z = exp(x);
79    J w = log(z);
80    VL << "z = " << z;
81    VL << "w = " << w;
82    ExpectJetsClose(w, x);
83  }
84
85  { // Check that (x * y) / x == y.
86    J z = x * y;
87    J w = z / x;
88    VL << "z = " << z;
89    VL << "w = " << w;
90    ExpectJetsClose(w, y);
91  }
92
93  { // Check that sqrt(x * x) == x.
94    J z = x * x;
95    J w = sqrt(z);
96    VL << "z = " << z;
97    VL << "w = " << w;
98    ExpectJetsClose(w, x);
99  }
100
101  { // Check that sqrt(y) * sqrt(y) == y.
102    J z = sqrt(y);
103    J w = z * z;
104    VL << "z = " << z;
105    VL << "w = " << w;
106    ExpectJetsClose(w, y);
107  }
108
109  { // Check that cos(2*x) = cos(x)^2 - sin(x)^2
110    J z = cos(J(2.0) * x);
111    J w = cos(x)*cos(x) - sin(x)*sin(x);
112    VL << "z = " << z;
113    VL << "w = " << w;
114    ExpectJetsClose(w, z);
115  }
116
117  { // Check that sin(2*x) = 2*cos(x)*sin(x)
118    J z = sin(J(2.0) * x);
119    J w = J(2.0)*cos(x)*sin(x);
120    VL << "z = " << z;
121    VL << "w = " << w;
122    ExpectJetsClose(w, z);
123  }
124
125  { // Check that cos(x)*cos(x) + sin(x)*sin(x) = 1
126    J z = cos(x) * cos(x);
127    J w = sin(x) * sin(x);
128    VL << "z = " << z;
129    VL << "w = " << w;
130    ExpectJetsClose(z + w, J(1.0));
131  }
132
133  { // Check that atan2(r*sin(t), r*cos(t)) = t.
134    J t = MakeJet(0.7, -0.3, +1.5);
135    J r = MakeJet(2.3, 0.13, -2.4);
136    VL << "t = " << t;
137    VL << "r = " << r;
138
139    J u = atan2(r * sin(t), r * cos(t));
140    VL << "u = " << u;
141
142    ExpectJetsClose(u, t);
143  }
144
145  { // Check that tan(x) = sin(x) / cos(x).
146    J z = tan(x);
147    J w = sin(x) / cos(x);
148    VL << "z = " << z;
149    VL << "w = " << w;
150    ExpectJetsClose(z, w);
151  }
152
153  { // Check that tan(atan(x)) = x.
154    J z = tan(atan(x));
155    J w = x;
156    VL << "z = " << z;
157    VL << "w = " << w;
158    ExpectJetsClose(z, w);
159  }
160
161  { // Check that cosh(x)*cosh(x) - sinh(x)*sinh(x) = 1
162    J z = cosh(x) * cosh(x);
163    J w = sinh(x) * sinh(x);
164    VL << "z = " << z;
165    VL << "w = " << w;
166    ExpectJetsClose(z - w, J(1.0));
167  }
168
169  { // Check that tanh(x + y) = (tanh(x) + tanh(y)) / (1 + tanh(x) tanh(y))
170    J z = tanh(x + y);
171    J w = (tanh(x) + tanh(y)) / (J(1.0) + tanh(x) * tanh(y));
172    VL << "z = " << z;
173    VL << "w = " << w;
174    ExpectJetsClose(z, w);
175  }
176
177  { // Check that pow(x, 1) == x.
178    VL << "x = " << x;
179
180    J u = pow(x, 1.);
181    VL << "u = " << u;
182
183    ExpectJetsClose(x, u);
184  }
185
186  { // Check that pow(x, 1) == x.
187    J y = MakeJet(1, 0.0, 0.0);
188    VL << "x = " << x;
189    VL << "y = " << y;
190
191    J u = pow(x, y);
192    VL << "u = " << u;
193
194    ExpectJetsClose(x, u);
195  }
196
197  { // Check that pow(e, log(x)) == x.
198    J logx = log(x);
199
200    VL << "x = " << x;
201    VL << "y = " << y;
202
203    J u = pow(kE, logx);
204    VL << "u = " << u;
205
206    ExpectJetsClose(x, u);
207  }
208
209  { // Check that pow(e, log(x)) == x.
210    J logx = log(x);
211    J e = MakeJet(kE, 0., 0.);
212    VL << "x = " << x;
213    VL << "log(x) = " << logx;
214
215    J u = pow(e, logx);
216    VL << "u = " << u;
217
218    ExpectJetsClose(x, u);
219  }
220
221  { // Check that pow(e, log(x)) == x.
222    J logx = log(x);
223    J e = MakeJet(kE, 0., 0.);
224    VL << "x = " << x;
225    VL << "logx = " << logx;
226
227    J u = pow(e, logx);
228    VL << "u = " << u;
229
230    ExpectJetsClose(x, u);
231  }
232
233  { // Check that pow(x,y) = exp(y*log(x)).
234    J logx = log(x);
235    J e = MakeJet(kE, 0., 0.);
236    VL << "x = " << x;
237    VL << "logx = " << logx;
238
239    J u = pow(e, y*logx);
240    J v = pow(x, y);
241    VL << "u = " << u;
242    VL << "v = " << v;
243
244    ExpectJetsClose(v, u);
245  }
246
247  { // Check that 1 + x == x + 1.
248    J a = x + 1.0;
249    J b = 1.0 + x;
250
251    ExpectJetsClose(a, b);
252  }
253
254  { // Check that 1 - x == -(x - 1).
255    J a = 1.0 - x;
256    J b = -(x - 1.0);
257
258    ExpectJetsClose(a, b);
259  }
260
261  { // Check that x/s == x*s.
262    J a = x / 5.0;
263    J b = x * 5.0;
264
265    ExpectJetsClose(5.0 * a, b / 5.0);
266  }
267
268  { // Check that x / y == 1 / (y / x).
269    J a = x / y;
270    J b = 1.0 / (y / x);
271    VL << "a = " << a;
272    VL << "b = " << b;
273
274    ExpectJetsClose(a, b);
275  }
276
277  { // Check that abs(-x * x) == sqrt(x * x).
278    ExpectJetsClose(abs(-x), sqrt(x * x));
279  }
280
281  { // Check that cos(acos(x)) == x.
282    J a = MakeJet(0.1, -2.7, 1e-3);
283    ExpectJetsClose(cos(acos(a)), a);
284    ExpectJetsClose(acos(cos(a)), a);
285
286    J b = MakeJet(0.6,  0.5, 1e+2);
287    ExpectJetsClose(cos(acos(b)), b);
288    ExpectJetsClose(acos(cos(b)), b);
289  }
290
291  { // Check that sin(asin(x)) == x.
292    J a = MakeJet(0.1, -2.7, 1e-3);
293    ExpectJetsClose(sin(asin(a)), a);
294    ExpectJetsClose(asin(sin(a)), a);
295
296    J b = MakeJet(0.4,  0.5, 1e+2);
297    ExpectJetsClose(sin(asin(b)), b);
298    ExpectJetsClose(asin(sin(b)), b);
299  }
300}
301
302TEST(Jet, JetsInEigenMatrices) {
303  J x = MakeJet(2.3, -2.7, 1e-3);
304  J y = MakeJet(1.7,  0.5, 1e+2);
305  J z = MakeJet(5.3, -4.7, 1e-3);
306  J w = MakeJet(9.7,  1.5, 10.1);
307
308  Eigen::Matrix<J, 2, 2> M;
309  Eigen::Matrix<J, 2, 1> v, r1, r2;
310
311  M << x, y, z, w;
312  v << x, z;
313
314  // Check that M * v == (v^T * M^T)^T
315  r1 = M * v;
316  r2 = (v.transpose() * M.transpose()).transpose();
317
318  ExpectJetsClose(r1(0), r2(0));
319  ExpectJetsClose(r1(1), r2(1));
320}
321
322TEST(JetTraitsTest, ClassificationMixed) {
323  Jet<double, 3> a(5.5, 0);
324  a.v[0] = std::numeric_limits<double>::quiet_NaN();
325  a.v[1] = std::numeric_limits<double>::infinity();
326  a.v[2] = -std::numeric_limits<double>::infinity();
327  EXPECT_FALSE(IsFinite(a));
328  EXPECT_FALSE(IsNormal(a));
329  EXPECT_TRUE(IsInfinite(a));
330  EXPECT_TRUE(IsNaN(a));
331}
332
333TEST(JetTraitsTest, ClassificationNaN) {
334  Jet<double, 3> a(5.5, 0);
335  a.v[0] = std::numeric_limits<double>::quiet_NaN();
336  a.v[1] = 0.0;
337  a.v[2] = 0.0;
338  EXPECT_FALSE(IsFinite(a));
339  EXPECT_FALSE(IsNormal(a));
340  EXPECT_FALSE(IsInfinite(a));
341  EXPECT_TRUE(IsNaN(a));
342}
343
344TEST(JetTraitsTest, ClassificationInf) {
345  Jet<double, 3> a(5.5, 0);
346  a.v[0] = std::numeric_limits<double>::infinity();
347  a.v[1] = 0.0;
348  a.v[2] = 0.0;
349  EXPECT_FALSE(IsFinite(a));
350  EXPECT_FALSE(IsNormal(a));
351  EXPECT_TRUE(IsInfinite(a));
352  EXPECT_FALSE(IsNaN(a));
353}
354
355TEST(JetTraitsTest, ClassificationFinite) {
356  Jet<double, 3> a(5.5, 0);
357  a.v[0] = 100.0;
358  a.v[1] = 1.0;
359  a.v[2] = 3.14159;
360  EXPECT_TRUE(IsFinite(a));
361  EXPECT_TRUE(IsNormal(a));
362  EXPECT_FALSE(IsInfinite(a));
363  EXPECT_FALSE(IsNaN(a));
364}
365
366}  // namespace internal
367}  // namespace ceres
368