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
44 namespace ceres {
45 namespace internal {
46
47 const double kE = 2.71828182845904523536;
48
49 typedef Jet<double, 2> J;
50
51 // Convenient shorthand for making a jet.
MakeJet(double a,double v0,double v1)52 J 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.
61 double const kTolerance = 1e-13;
62
ExpectJetsClose(const J & x,const J & y)63 void 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
TEST(Jet,Jet)69 TEST(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
TEST(Jet,JetsInEigenMatrices)302 TEST(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
TEST(JetTraitsTest,ClassificationMixed)322 TEST(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
TEST(JetTraitsTest,ClassificationNaN)333 TEST(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
TEST(JetTraitsTest,ClassificationInf)344 TEST(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
TEST(JetTraitsTest,ClassificationFinite)355 TEST(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