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 (at) 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. 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 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 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 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 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 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 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 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
