internal/ceres/jet_test.cc

   1 // Ceres Solver - A fast non-linear least squares minimizer
   2 // Copyright 2015 Google Inc. All rights reserved.
   3 // http://ceres-solver.org/
   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.
  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 pow(0, y) == 0 for y > 1, with both arguments Jets.
 248     // This tests special case handling inside pow().
 249     J a = MakeJet(0, 1, 2);
 250     J b = MakeJet(2, 3, 4);
 251     VL << "a = " << a;
 252     VL << "b = " << b;
 253
 254     J c = pow(a, b);
 255     VL << "a^b = " << c;
 256     ExpectJetsClose(c, MakeJet(0, 0, 0));
 257   }
 258
 259   { // Check that pow(0, y) == 0 for y == 1, with both arguments Jets.
 260     // This tests special case handling inside pow().
 261     J a = MakeJet(0, 1, 2);
 262     J b = MakeJet(1, 3, 4);
 263     VL << "a = " << a;
 264     VL << "b = " << b;
 265
 266     J c = pow(a, b);
 267     VL << "a^b = " << c;
 268     ExpectJetsClose(c, MakeJet(0, 1, 2));
 269   }
 270
 271   { // Check that pow(0, <1) is not finite, with both arguments Jets.
 272     for (int i = 1; i < 10; i++) {
 273       J a = MakeJet(0, 1, 2);
 274       J b = MakeJet(i*0.1, 3, 4);       // b = 0.1 ... 0.9
 275       VL << "a = " << a;
 276       VL << "b = " << b;
 277
 278       J c = pow(a, b);
 279       VL << "a^b = " << c;
 280       EXPECT_EQ(c.a, 0.0);
 281       EXPECT_FALSE(IsFinite(c.v[0]));
 282       EXPECT_FALSE(IsFinite(c.v[1]));
 283     }
 284     for (int i = -10; i < 0; i++) {
 285       J a = MakeJet(0, 1, 2);
 286       J b = MakeJet(i*0.1, 3, 4);       // b = -1,-0.9 ... -0.1
 287       VL << "a = " << a;
 288       VL << "b = " << b;
 289
 290       J c = pow(a, b);
 291       VL << "a^b = " << c;
 292       EXPECT_FALSE(IsFinite(c.a));
 293       EXPECT_FALSE(IsFinite(c.v[0]));
 294       EXPECT_FALSE(IsFinite(c.v[1]));
 295     }
 296
 297     {
 298       // The special case of 0^0 = 1 defined by the C standard.
 299       J a = MakeJet(0, 1, 2);
 300       J b = MakeJet(0, 3, 4);
 301       VL << "a = " << a;
 302       VL << "b = " << b;
 303
 304       J c = pow(a, b);
 305       VL << "a^b = " << c;
 306       EXPECT_EQ(c.a, 1.0);
 307       EXPECT_FALSE(IsFinite(c.v[0]));
 308       EXPECT_FALSE(IsFinite(c.v[1]));
 309     }
 310   }
 311
 312   { // Check that pow(<0, b) is correct for integer b.
 313     // This tests special case handling inside pow().
 314     J a = MakeJet(-1.5, 3, 4);
 315
 316     // b integer:
 317     for (int i = -10; i <= 10; i++) {
 318       J b = MakeJet(i, 0, 5);
 319       VL << "a = " << a;
 320       VL << "b = " << b;
 321
 322       J c = pow(a, b);
 323       VL << "a^b = " << c;
 324       ExpectClose(c.a, pow(-1.5, i), kTolerance);
 325       EXPECT_TRUE(IsFinite(c.v[0]));
 326       EXPECT_FALSE(IsFinite(c.v[1]));
 327       ExpectClose(c.v[0], i * pow(-1.5, i - 1) * 3.0, kTolerance);
 328     }
 329   }
 330
 331   { // Check that pow(<0, b) is correct for noninteger b.
 332     // This tests special case handling inside pow().
 333     J a = MakeJet(-1.5, 3, 4);
 334     J b = MakeJet(-2.5, 0, 5);
 335     VL << "a = " << a;
 336     VL << "b = " << b;
 337
 338     J c = pow(a, b);
 339     VL << "a^b = " << c;
 340     EXPECT_FALSE(IsFinite(c.a));
 341     EXPECT_FALSE(IsFinite(c.v[0]));
 342     EXPECT_FALSE(IsFinite(c.v[1]));
 343   }
 344
 345   {
 346     // Check that pow(0,y) == 0 for y == 2, with the second argument a
 347     // Jet.  This tests special case handling inside pow().
 348     double a = 0;
 349     J b = MakeJet(2, 3, 4);
 350     VL << "a = " << a;
 351     VL << "b = " << b;
 352
 353     J c = pow(a, b);
 354     VL << "a^b = " << c;
 355     ExpectJetsClose(c, MakeJet(0, 0, 0));
 356   }
 357
 358   {
 359     // Check that pow(<0,y) is correct for integer y. This tests special case
 360     // handling inside pow().
 361     double a = -1.5;
 362     for (int i = -10; i <= 10; i++) {
 363       J b = MakeJet(i, 3, 0);
 364       VL << "a = " << a;
 365       VL << "b = " << b;
 366
 367       J c = pow(a, b);
 368       VL << "a^b = " << c;
 369       ExpectClose(c.a, pow(-1.5, i), kTolerance);
 370       EXPECT_FALSE(IsFinite(c.v[0]));
 371       EXPECT_TRUE(IsFinite(c.v[1]));
 372       ExpectClose(c.v[1], 0, kTolerance);
 373     }
 374   }
 375
 376   {
 377     // Check that pow(<0,y) is correct for noninteger y. This tests special
 378     // case handling inside pow().
 379     double a = -1.5;
 380     J b = MakeJet(-3.14, 3, 0);
 381     VL << "a = " << a;
 382     VL << "b = " << b;
 383
 384     J c = pow(a, b);
 385     VL << "a^b = " << c;
 386     EXPECT_FALSE(IsFinite(c.a));
 387     EXPECT_FALSE(IsFinite(c.v[0]));
 388     EXPECT_FALSE(IsFinite(c.v[1]));
 389   }
 390
 391   { // Check that 1 + x == x + 1.
 392     J a = x + 1.0;
 393     J b = 1.0 + x;
 394     J c = x;
 395     c += 1.0;
 396
 397     ExpectJetsClose(a, b);
 398     ExpectJetsClose(a, c);
 399   }
 400
 401   { // Check that 1 - x == -(x - 1).
 402     J a = 1.0 - x;
 403     J b = -(x - 1.0);
 404     J c = x;
 405     c -= 1.0;
 406
 407     ExpectJetsClose(a, b);
 408     ExpectJetsClose(a, -c);
 409   }
 410
 411   { // Check that (x/s)*s == (x*s)/s.
 412     J a = x / 5.0;
 413     J b = x * 5.0;
 414     J c = x;
 415     c /= 5.0;
 416     J d = x;
 417     d *= 5.0;
 418
 419     ExpectJetsClose(5.0 * a, b / 5.0);
 420     ExpectJetsClose(a, c);
 421     ExpectJetsClose(b, d);
 422   }
 423
 424   { // Check that x / y == 1 / (y / x).
 425     J a = x / y;
 426     J b = 1.0 / (y / x);
 427     VL << "a = " << a;
 428     VL << "b = " << b;
 429
 430     ExpectJetsClose(a, b);
 431   }
 432
 433   { // Check that abs(-x * x) == sqrt(x * x).
 434     ExpectJetsClose(abs(-x), sqrt(x * x));
 435   }
 436
 437   { // Check that cos(acos(x)) == x.
 438     J a = MakeJet(0.1, -2.7, 1e-3);
 439     ExpectJetsClose(cos(acos(a)), a);
 440     ExpectJetsClose(acos(cos(a)), a);
 441
 442     J b = MakeJet(0.6,  0.5, 1e+2);
 443     ExpectJetsClose(cos(acos(b)), b);
 444     ExpectJetsClose(acos(cos(b)), b);
 445   }
 446
 447   { // Check that sin(asin(x)) == x.
 448     J a = MakeJet(0.1, -2.7, 1e-3);
 449     ExpectJetsClose(sin(asin(a)), a);
 450     ExpectJetsClose(asin(sin(a)), a);
 451
 452     J b = MakeJet(0.4,  0.5, 1e+2);
 453     ExpectJetsClose(sin(asin(b)), b);
 454     ExpectJetsClose(asin(sin(b)), b);
 455   }
 456
 457   {
 458     J zero = J(0.0);
 459
 460     // Check that J0(0) == 1.
 461     ExpectJetsClose(BesselJ0(zero), J(1.0));
 462
 463     // Check that J1(0) == 0.
 464     ExpectJetsClose(BesselJ1(zero), zero);
 465
 466     // Check that J2(0) == 0.
 467     ExpectJetsClose(BesselJn(2, zero), zero);
 468
 469     // Check that J3(0) == 0.
 470     ExpectJetsClose(BesselJn(3, zero), zero);
 471
 472     J z = MakeJet(0.1, -2.7, 1e-3);
 473
 474     // Check that J0(z) == Jn(0,z).
 475     ExpectJetsClose(BesselJ0(z), BesselJn(0, z));
 476
 477     // Check that J1(z) == Jn(1,z).
 478     ExpectJetsClose(BesselJ1(z), BesselJn(1, z));
 479
 480     // Check that J0(z)+J2(z) == (2/z)*J1(z).
 481     // See formula http://dlmf.nist.gov/10.6.E1
 482     ExpectJetsClose(BesselJ0(z) + BesselJn(2, z), (2.0 / z) * BesselJ1(z));
 483   }
 484
 485   { // Check that floor of a positive number works.
 486     J a = MakeJet(0.1, -2.7, 1e-3);
 487     J b = floor(a);
 488     J expected = MakeJet(floor(a.a), 0.0, 0.0);
 489     EXPECT_EQ(expected, b);
 490   }
 491
 492   { // Check that floor of a negative number works.
 493     J a = MakeJet(-1.1, -2.7, 1e-3);
 494     J b = floor(a);
 495     J expected = MakeJet(floor(a.a), 0.0, 0.0);
 496     EXPECT_EQ(expected, b);
 497   }
 498
 499   { // Check that floor of a positive number works.
 500     J a = MakeJet(10.123, -2.7, 1e-3);
 501     J b = floor(a);
 502     J expected = MakeJet(floor(a.a), 0.0, 0.0);
 503     EXPECT_EQ(expected, b);
 504   }
 505
 506   { // Check that ceil of a positive number works.
 507     J a = MakeJet(0.1, -2.7, 1e-3);
 508     J b = ceil(a);
 509     J expected = MakeJet(ceil(a.a), 0.0, 0.0);
 510     EXPECT_EQ(expected, b);
 511   }
 512
 513   { // Check that ceil of a negative number works.
 514     J a = MakeJet(-1.1, -2.7, 1e-3);
 515     J b = ceil(a);
 516     J expected = MakeJet(ceil(a.a), 0.0, 0.0);
 517     EXPECT_EQ(expected, b);
 518   }
 519
 520   { // Check that ceil of a positive number works.
 521     J a = MakeJet(10.123, -2.7, 1e-3);
 522     J b = ceil(a);
 523     J expected = MakeJet(ceil(a.a), 0.0, 0.0);
 524     EXPECT_EQ(expected, b);
 525   }
 526 }
 527
 528 TEST(Jet, JetsInEigenMatrices) {
 529   J x = MakeJet(2.3, -2.7, 1e-3);
 530   J y = MakeJet(1.7,  0.5, 1e+2);
 531   J z = MakeJet(5.3, -4.7, 1e-3);
 532   J w = MakeJet(9.7,  1.5, 10.1);
 533
 534   Eigen::Matrix<J, 2, 2> M;
 535   Eigen::Matrix<J, 2, 1> v, r1, r2;
 536
 537   M << x, y, z, w;
 538   v << x, z;
 539
 540   // Check that M * v == (v^T * M^T)^T
 541   r1 = M * v;
 542   r2 = (v.transpose() * M.transpose()).transpose();
 543
 544   ExpectJetsClose(r1(0), r2(0));
 545   ExpectJetsClose(r1(1), r2(1));
 546 }
 547
 548 TEST(JetTraitsTest, ClassificationMixed) {
 549   Jet<double, 3> a(5.5, 0);
 550   a.v[0] = std::numeric_limits<double>::quiet_NaN();
 551   a.v[1] = std::numeric_limits<double>::infinity();
 552   a.v[2] = -std::numeric_limits<double>::infinity();
 553   EXPECT_FALSE(IsFinite(a));
 554   EXPECT_FALSE(IsNormal(a));
 555   EXPECT_TRUE(IsInfinite(a));
 556   EXPECT_TRUE(IsNaN(a));
 557 }
 558
 559 TEST(JetTraitsTest, ClassificationNaN) {
 560   Jet<double, 3> a(5.5, 0);
 561   a.v[0] = std::numeric_limits<double>::quiet_NaN();
 562   a.v[1] = 0.0;
 563   a.v[2] = 0.0;
 564   EXPECT_FALSE(IsFinite(a));
 565   EXPECT_FALSE(IsNormal(a));
 566   EXPECT_FALSE(IsInfinite(a));
 567   EXPECT_TRUE(IsNaN(a));
 568 }
 569
 570 TEST(JetTraitsTest, ClassificationInf) {
 571   Jet<double, 3> a(5.5, 0);
 572   a.v[0] = std::numeric_limits<double>::infinity();
 573   a.v[1] = 0.0;
 574   a.v[2] = 0.0;
 575   EXPECT_FALSE(IsFinite(a));
 576   EXPECT_FALSE(IsNormal(a));
 577   EXPECT_TRUE(IsInfinite(a));
 578   EXPECT_FALSE(IsNaN(a));
 579 }
 580
 581 TEST(JetTraitsTest, ClassificationFinite) {
 582   Jet<double, 3> a(5.5, 0);
 583   a.v[0] = 100.0;
 584   a.v[1] = 1.0;
 585   a.v[2] = 3.14159;
 586   EXPECT_TRUE(IsFinite(a));
 587   EXPECT_TRUE(IsNormal(a));
 588   EXPECT_FALSE(IsInfinite(a));
 589   EXPECT_FALSE(IsNaN(a));
 590 }
 591
 592 // ScalarBinaryOpTraits is only supported on Eigen versions >= 3.3
 593 #if EIGEN_VERSION_AT_LEAST(3, 3, 0)
 594 TEST(JetTraitsTest, MatrixScalarUnaryOps) {
 595   const J x = MakeJet(2.3, -2.7, 1e-3);
 596   const J y = MakeJet(1.7,  0.5, 1e+2);
 597   Eigen::Matrix<J, 2, 1> a;
 598   a << x, y;
 599
 600   const J sum = a.sum();
 601   const J sum2 = a(0) + a(1);
 602   ExpectJetsClose(sum, sum2);
 603 }
 604
 605 TEST(JetTraitsTest, MatrixScalarBinaryOps) {
 606   const J x = MakeJet(2.3, -2.7, 1e-3);
 607   const J y = MakeJet(1.7,  0.5, 1e+2);
 608   const J z = MakeJet(5.3, -4.7, 1e-3);
 609   const J w = MakeJet(9.7,  1.5, 10.1);
 610
 611   Eigen::Matrix<J, 2, 2> M;
 612   Eigen::Vector2d v;
 613
 614   M << x, y, z, w;
 615   v << 0.6, -2.1;
 616
 617   // Check that M * v == M * v.cast<J>().
 618   const Eigen::Matrix<J, 2, 1> r1 = M * v;
 619   const Eigen::Matrix<J, 2, 1> r2 = M * v.cast<J>();
 620
 621   ExpectJetsClose(r1(0), r2(0));
 622   ExpectJetsClose(r1(1), r2(1));
 623
 624   // Check that M * a == M * T(a).
 625   const double a = 3.1;
 626   const Eigen::Matrix<J, 2, 2> r3 = M * a;
 627   const Eigen::Matrix<J, 2, 2> r4 = M * J(a);
 628
 629   ExpectJetsClose(r3(0, 0), r4(0, 0));
 630   ExpectJetsClose(r3(1, 0), r4(1, 0));
 631   ExpectJetsClose(r3(0, 1), r4(0, 1));
 632   ExpectJetsClose(r3(1, 1), r4(1, 1));
 633 }
 634
 635 TEST(JetTraitsTest, ArrayScalarUnaryOps) {
 636   const J x = MakeJet(2.3, -2.7, 1e-3);
 637   const J y = MakeJet(1.7,  0.5, 1e+2);
 638   Eigen::Array<J, 2, 1> a;
 639   a << x, y;
 640
 641   const J sum = a.sum();
 642   const J sum2 = a(0) + a(1);
 643   ExpectJetsClose(sum, sum2);
 644 }
 645
 646 TEST(JetTraitsTest, ArrayScalarBinaryOps) {
 647   const J x = MakeJet(2.3, -2.7, 1e-3);
 648   const J y = MakeJet(1.7,  0.5, 1e+2);
 649
 650   Eigen::Array<J, 2, 1> a;
 651   Eigen::Array2d b;
 652
 653   a << x, y;
 654   b << 0.6, -2.1;
 655
 656   // Check that a * b == a * b.cast<T>()
 657   const Eigen::Array<J, 2, 1> r1 = a * b;
 658   const Eigen::Array<J, 2, 1> r2 = a * b.cast<J>();
 659
 660   ExpectJetsClose(r1(0), r2(0));
 661   ExpectJetsClose(r1(1), r2(1));
 662
 663   // Check that a * c == a * T(c).
 664   const double c = 3.1;
 665   const Eigen::Array<J, 2, 1> r3 = a * c;
 666   const Eigen::Array<J, 2, 1> r4 = a * J(c);
 667
 668   ExpectJetsClose(r3(0), r3(0));
 669   ExpectJetsClose(r4(1), r4(1));
 670 }
 671 #endif   // EIGEN_VERSION_AT_LEAST(3, 3, 0)
 672
 673 }  // namespace internal
 674 }  // namespace ceres