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// A simple implementation of N-dimensional dual numbers, for automatically 32// computing exact derivatives of functions. 33// 34// While a complete treatment of the mechanics of automatic differentation is 35// beyond the scope of this header (see 36// http://en.wikipedia.org/wiki/Automatic_differentiation for details), the 37// basic idea is to extend normal arithmetic with an extra element, "e," often 38// denoted with the greek symbol epsilon, such that e != 0 but e^2 = 0. Dual 39// numbers are extensions of the real numbers analogous to complex numbers: 40// whereas complex numbers augment the reals by introducing an imaginary unit i 41// such that i^2 = -1, dual numbers introduce an "infinitesimal" unit e such 42// that e^2 = 0. Dual numbers have two components: the "real" component and the 43// "infinitesimal" component, generally written as x + y*e. Surprisingly, this 44// leads to a convenient method for computing exact derivatives without needing 45// to manipulate complicated symbolic expressions. 46// 47// For example, consider the function 48// 49// f(x) = x^2 , 50// 51// evaluated at 10. Using normal arithmetic, f(10) = 100, and df/dx(10) = 20. 52// Next, augument 10 with an infinitesimal to get: 53// 54// f(10 + e) = (10 + e)^2 55// = 100 + 2 * 10 * e + e^2 56// = 100 + 20 * e -+- 57// -- | 58// | +--- This is zero, since e^2 = 0 59// | 60// +----------------- This is df/dx! 61// 62// Note that the derivative of f with respect to x is simply the infinitesimal 63// component of the value of f(x + e). So, in order to take the derivative of 64// any function, it is only necessary to replace the numeric "object" used in 65// the function with one extended with infinitesimals. The class Jet, defined in 66// this header, is one such example of this, where substitution is done with 67// templates. 68// 69// To handle derivatives of functions taking multiple arguments, different 70// infinitesimals are used, one for each variable to take the derivative of. For 71// example, consider a scalar function of two scalar parameters x and y: 72// 73// f(x, y) = x^2 + x * y 74// 75// Following the technique above, to compute the derivatives df/dx and df/dy for 76// f(1, 3) involves doing two evaluations of f, the first time replacing x with 77// x + e, the second time replacing y with y + e. 78// 79// For df/dx: 80// 81// f(1 + e, y) = (1 + e)^2 + (1 + e) * 3 82// = 1 + 2 * e + 3 + 3 * e 83// = 4 + 5 * e 84// 85// --> df/dx = 5 86// 87// For df/dy: 88// 89// f(1, 3 + e) = 1^2 + 1 * (3 + e) 90// = 1 + 3 + e 91// = 4 + e 92// 93// --> df/dy = 1 94// 95// To take the gradient of f with the implementation of dual numbers ("jets") in 96// this file, it is necessary to create a single jet type which has components 97// for the derivative in x and y, and passing them to a templated version of f: 98// 99// template<typename T> 100// T f(const T &x, const T &y) { 101// return x * x + x * y; 102// } 103// 104// // The "2" means there should be 2 dual number components. 105// Jet<double, 2> x(0); // Pick the 0th dual number for x. 106// Jet<double, 2> y(1); // Pick the 1st dual number for y. 107// Jet<double, 2> z = f(x, y); 108// 109// LG << "df/dx = " << z.a[0] 110// << "df/dy = " << z.a[1]; 111// 112// Most users should not use Jet objects directly; a wrapper around Jet objects, 113// which makes computing the derivative, gradient, or jacobian of templated 114// functors simple, is in autodiff.h. Even autodiff.h should not be used 115// directly; instead autodiff_cost_function.h is typically the file of interest. 116// 117// For the more mathematically inclined, this file implements first-order 118// "jets". A 1st order jet is an element of the ring 119// 120// T[N] = T[t_1, ..., t_N] / (t_1, ..., t_N)^2 121// 122// which essentially means that each jet consists of a "scalar" value 'a' from T 123// and a 1st order perturbation vector 'v' of length N: 124// 125// x = a + \sum_i v[i] t_i 126// 127// A shorthand is to write an element as x = a + u, where u is the pertubation. 128// Then, the main point about the arithmetic of jets is that the product of 129// perturbations is zero: 130// 131// (a + u) * (b + v) = ab + av + bu + uv 132// = ab + (av + bu) + 0 133// 134// which is what operator* implements below. Addition is simpler: 135// 136// (a + u) + (b + v) = (a + b) + (u + v). 137// 138// The only remaining question is how to evaluate the function of a jet, for 139// which we use the chain rule: 140// 141// f(a + u) = f(a) + f'(a) u 142// 143// where f'(a) is the (scalar) derivative of f at a. 144// 145// By pushing these things through sufficiently and suitably templated 146// functions, we can do automatic differentiation. Just be sure to turn on 147// function inlining and common-subexpression elimination, or it will be very 148// slow! 149// 150// WARNING: Most Ceres users should not directly include this file or know the 151// details of how jets work. Instead the suggested method for automatic 152// derivatives is to use autodiff_cost_function.h, which is a wrapper around 153// both jets.h and autodiff.h to make taking derivatives of cost functions for 154// use in Ceres easier. 155 156#ifndef CERES_PUBLIC_JET_H_ 157#define CERES_PUBLIC_JET_H_ 158 159#include <cmath> 160#include <iosfwd> 161#include <iostream> // NOLINT 162#include <string> 163 164#include "Eigen/Core" 165#include "ceres/fpclassify.h" 166 167namespace ceres { 168 169template <typename T, int N> 170struct Jet { 171 enum { DIMENSION = N }; 172 173 // Default-construct "a" because otherwise this can lead to false errors about 174 // uninitialized uses when other classes relying on default constructed T 175 // (where T is a Jet<T, N>). This usually only happens in opt mode. Note that 176 // the C++ standard mandates that e.g. default constructed doubles are 177 // initialized to 0.0; see sections 8.5 of the C++03 standard. 178 Jet() : a() { 179 v.setZero(); 180 } 181 182 // Constructor from scalar: a + 0. 183 explicit Jet(const T& value) { 184 a = value; 185 v.setZero(); 186 } 187 188 // Constructor from scalar plus variable: a + t_i. 189 Jet(const T& value, int k) { 190 a = value; 191 v.setZero(); 192 v[k] = T(1.0); 193 } 194 195 // Compound operators 196 Jet<T, N>& operator+=(const Jet<T, N> &y) { 197 *this = *this + y; 198 return *this; 199 } 200 201 Jet<T, N>& operator-=(const Jet<T, N> &y) { 202 *this = *this - y; 203 return *this; 204 } 205 206 Jet<T, N>& operator*=(const Jet<T, N> &y) { 207 *this = *this * y; 208 return *this; 209 } 210 211 Jet<T, N>& operator/=(const Jet<T, N> &y) { 212 *this = *this / y; 213 return *this; 214 } 215 216 // The scalar part. 217 T a; 218 219 // The infinitesimal part. 220 // 221 // Note the Eigen::DontAlign bit is needed here because this object 222 // gets allocated on the stack and as part of other arrays and 223 // structs. Forcing the right alignment there is the source of much 224 // pain and suffering. Even if that works, passing Jets around to 225 // functions by value has problems because the C++ ABI does not 226 // guarantee alignment for function arguments. 227 // 228 // Setting the DontAlign bit prevents Eigen from using SSE for the 229 // various operations on Jets. This is a small performance penalty 230 // since the AutoDiff code will still expose much of the code as 231 // statically sized loops to the compiler. But given the subtle 232 // issues that arise due to alignment, especially when dealing with 233 // multiple platforms, it seems to be a trade off worth making. 234 Eigen::Matrix<T, N, 1, Eigen::DontAlign> v; 235}; 236 237// Unary + 238template<typename T, int N> inline 239Jet<T, N> const& operator+(const Jet<T, N>& f) { 240 return f; 241} 242 243// TODO(keir): Try adding __attribute__((always_inline)) to these functions to 244// see if it causes a performance increase. 245 246// Unary - 247template<typename T, int N> inline 248Jet<T, N> operator-(const Jet<T, N>&f) { 249 Jet<T, N> g; 250 g.a = -f.a; 251 g.v = -f.v; 252 return g; 253} 254 255// Binary + 256template<typename T, int N> inline 257Jet<T, N> operator+(const Jet<T, N>& f, 258 const Jet<T, N>& g) { 259 Jet<T, N> h; 260 h.a = f.a + g.a; 261 h.v = f.v + g.v; 262 return h; 263} 264 265// Binary + with a scalar: x + s 266template<typename T, int N> inline 267Jet<T, N> operator+(const Jet<T, N>& f, T s) { 268 Jet<T, N> h; 269 h.a = f.a + s; 270 h.v = f.v; 271 return h; 272} 273 274// Binary + with a scalar: s + x 275template<typename T, int N> inline 276Jet<T, N> operator+(T s, const Jet<T, N>& f) { 277 Jet<T, N> h; 278 h.a = f.a + s; 279 h.v = f.v; 280 return h; 281} 282 283// Binary - 284template<typename T, int N> inline 285Jet<T, N> operator-(const Jet<T, N>& f, 286 const Jet<T, N>& g) { 287 Jet<T, N> h; 288 h.a = f.a - g.a; 289 h.v = f.v - g.v; 290 return h; 291} 292 293// Binary - with a scalar: x - s 294template<typename T, int N> inline 295Jet<T, N> operator-(const Jet<T, N>& f, T s) { 296 Jet<T, N> h; 297 h.a = f.a - s; 298 h.v = f.v; 299 return h; 300} 301 302// Binary - with a scalar: s - x 303template<typename T, int N> inline 304Jet<T, N> operator-(T s, const Jet<T, N>& f) { 305 Jet<T, N> h; 306 h.a = s - f.a; 307 h.v = -f.v; 308 return h; 309} 310 311// Binary * 312template<typename T, int N> inline 313Jet<T, N> operator*(const Jet<T, N>& f, 314 const Jet<T, N>& g) { 315 Jet<T, N> h; 316 h.a = f.a * g.a; 317 h.v = f.a * g.v + f.v * g.a; 318 return h; 319} 320 321// Binary * with a scalar: x * s 322template<typename T, int N> inline 323Jet<T, N> operator*(const Jet<T, N>& f, T s) { 324 Jet<T, N> h; 325 h.a = f.a * s; 326 h.v = f.v * s; 327 return h; 328} 329 330// Binary * with a scalar: s * x 331template<typename T, int N> inline 332Jet<T, N> operator*(T s, const Jet<T, N>& f) { 333 Jet<T, N> h; 334 h.a = f.a * s; 335 h.v = f.v * s; 336 return h; 337} 338 339// Binary / 340template<typename T, int N> inline 341Jet<T, N> operator/(const Jet<T, N>& f, 342 const Jet<T, N>& g) { 343 Jet<T, N> h; 344 // This uses: 345 // 346 // a + u (a + u)(b - v) (a + u)(b - v) 347 // ----- = -------------- = -------------- 348 // b + v (b + v)(b - v) b^2 349 // 350 // which holds because v*v = 0. 351 h.a = f.a / g.a; 352 h.v = (f.v - f.a / g.a * g.v) / g.a; 353 return h; 354} 355 356// Binary / with a scalar: s / x 357template<typename T, int N> inline 358Jet<T, N> operator/(T s, const Jet<T, N>& g) { 359 Jet<T, N> h; 360 h.a = s / g.a; 361 h.v = - s * g.v / (g.a * g.a); 362 return h; 363} 364 365// Binary / with a scalar: x / s 366template<typename T, int N> inline 367Jet<T, N> operator/(const Jet<T, N>& f, T s) { 368 Jet<T, N> h; 369 h.a = f.a / s; 370 h.v = f.v / s; 371 return h; 372} 373 374// Binary comparison operators for both scalars and jets. 375#define CERES_DEFINE_JET_COMPARISON_OPERATOR(op) \ 376template<typename T, int N> inline \ 377bool operator op(const Jet<T, N>& f, const Jet<T, N>& g) { \ 378 return f.a op g.a; \ 379} \ 380template<typename T, int N> inline \ 381bool operator op(const T& s, const Jet<T, N>& g) { \ 382 return s op g.a; \ 383} \ 384template<typename T, int N> inline \ 385bool operator op(const Jet<T, N>& f, const T& s) { \ 386 return f.a op s; \ 387} 388CERES_DEFINE_JET_COMPARISON_OPERATOR( < ) // NOLINT 389CERES_DEFINE_JET_COMPARISON_OPERATOR( <= ) // NOLINT 390CERES_DEFINE_JET_COMPARISON_OPERATOR( > ) // NOLINT 391CERES_DEFINE_JET_COMPARISON_OPERATOR( >= ) // NOLINT 392CERES_DEFINE_JET_COMPARISON_OPERATOR( == ) // NOLINT 393CERES_DEFINE_JET_COMPARISON_OPERATOR( != ) // NOLINT 394#undef CERES_DEFINE_JET_COMPARISON_OPERATOR 395 396// Pull some functions from namespace std. 397// 398// This is necessary because we want to use the same name (e.g. 'sqrt') for 399// double-valued and Jet-valued functions, but we are not allowed to put 400// Jet-valued functions inside namespace std. 401// 402// Missing: cosh, sinh, tanh, tan 403// TODO(keir): Switch to "using". 404inline double abs (double x) { return std::abs(x); } 405inline double log (double x) { return std::log(x); } 406inline double exp (double x) { return std::exp(x); } 407inline double sqrt (double x) { return std::sqrt(x); } 408inline double cos (double x) { return std::cos(x); } 409inline double acos (double x) { return std::acos(x); } 410inline double sin (double x) { return std::sin(x); } 411inline double asin (double x) { return std::asin(x); } 412inline double pow (double x, double y) { return std::pow(x, y); } 413inline double atan2(double y, double x) { return std::atan2(y, x); } 414 415// In general, f(a + h) ~= f(a) + f'(a) h, via the chain rule. 416 417// abs(x + h) ~= x + h or -(x + h) 418template <typename T, int N> inline 419Jet<T, N> abs(const Jet<T, N>& f) { 420 return f.a < T(0.0) ? -f : f; 421} 422 423// log(a + h) ~= log(a) + h / a 424template <typename T, int N> inline 425Jet<T, N> log(const Jet<T, N>& f) { 426 Jet<T, N> g; 427 g.a = log(f.a); 428 g.v = f.v / f.a; 429 return g; 430} 431 432// exp(a + h) ~= exp(a) + exp(a) h 433template <typename T, int N> inline 434Jet<T, N> exp(const Jet<T, N>& f) { 435 Jet<T, N> g; 436 g.a = exp(f.a); 437 g.v = g.a * f.v; 438 return g; 439} 440 441// sqrt(a + h) ~= sqrt(a) + h / (2 sqrt(a)) 442template <typename T, int N> inline 443Jet<T, N> sqrt(const Jet<T, N>& f) { 444 Jet<T, N> g; 445 g.a = sqrt(f.a); 446 g.v = f.v / (T(2.0) * g.a); 447 return g; 448} 449 450// cos(a + h) ~= cos(a) - sin(a) h 451template <typename T, int N> inline 452Jet<T, N> cos(const Jet<T, N>& f) { 453 Jet<T, N> g; 454 g.a = cos(f.a); 455 T sin_a = sin(f.a); 456 g.v = - sin_a * f.v; 457 return g; 458} 459 460// acos(a + h) ~= acos(a) - 1 / sqrt(1 - a^2) h 461template <typename T, int N> inline 462Jet<T, N> acos(const Jet<T, N>& f) { 463 Jet<T, N> g; 464 g.a = acos(f.a); 465 g.v = - T(1.0) / sqrt(T(1.0) - f.a * f.a) * f.v; 466 return g; 467} 468 469// sin(a + h) ~= sin(a) + cos(a) h 470template <typename T, int N> inline 471Jet<T, N> sin(const Jet<T, N>& f) { 472 Jet<T, N> g; 473 g.a = sin(f.a); 474 T cos_a = cos(f.a); 475 g.v = cos_a * f.v; 476 return g; 477} 478 479// asin(a + h) ~= asin(a) + 1 / sqrt(1 - a^2) h 480template <typename T, int N> inline 481Jet<T, N> asin(const Jet<T, N>& f) { 482 Jet<T, N> g; 483 g.a = asin(f.a); 484 g.v = T(1.0) / sqrt(T(1.0) - f.a * f.a) * f.v; 485 return g; 486} 487 488// Jet Classification. It is not clear what the appropriate semantics are for 489// these classifications. This picks that IsFinite and isnormal are "all" 490// operations, i.e. all elements of the jet must be finite for the jet itself 491// to be finite (or normal). For IsNaN and IsInfinite, the answer is less 492// clear. This takes a "any" approach for IsNaN and IsInfinite such that if any 493// part of a jet is nan or inf, then the entire jet is nan or inf. This leads 494// to strange situations like a jet can be both IsInfinite and IsNaN, but in 495// practice the "any" semantics are the most useful for e.g. checking that 496// derivatives are sane. 497 498// The jet is finite if all parts of the jet are finite. 499template <typename T, int N> inline 500bool IsFinite(const Jet<T, N>& f) { 501 if (!IsFinite(f.a)) { 502 return false; 503 } 504 for (int i = 0; i < N; ++i) { 505 if (!IsFinite(f.v[i])) { 506 return false; 507 } 508 } 509 return true; 510} 511 512// The jet is infinite if any part of the jet is infinite. 513template <typename T, int N> inline 514bool IsInfinite(const Jet<T, N>& f) { 515 if (IsInfinite(f.a)) { 516 return true; 517 } 518 for (int i = 0; i < N; i++) { 519 if (IsInfinite(f.v[i])) { 520 return true; 521 } 522 } 523 return false; 524} 525 526// The jet is NaN if any part of the jet is NaN. 527template <typename T, int N> inline 528bool IsNaN(const Jet<T, N>& f) { 529 if (IsNaN(f.a)) { 530 return true; 531 } 532 for (int i = 0; i < N; ++i) { 533 if (IsNaN(f.v[i])) { 534 return true; 535 } 536 } 537 return false; 538} 539 540// The jet is normal if all parts of the jet are normal. 541template <typename T, int N> inline 542bool IsNormal(const Jet<T, N>& f) { 543 if (!IsNormal(f.a)) { 544 return false; 545 } 546 for (int i = 0; i < N; ++i) { 547 if (!IsNormal(f.v[i])) { 548 return false; 549 } 550 } 551 return true; 552} 553 554// atan2(b + db, a + da) ~= atan2(b, a) + (- b da + a db) / (a^2 + b^2) 555// 556// In words: the rate of change of theta is 1/r times the rate of 557// change of (x, y) in the positive angular direction. 558template <typename T, int N> inline 559Jet<T, N> atan2(const Jet<T, N>& g, const Jet<T, N>& f) { 560 // Note order of arguments: 561 // 562 // f = a + da 563 // g = b + db 564 565 Jet<T, N> out; 566 567 out.a = atan2(g.a, f.a); 568 569 T const temp = T(1.0) / (f.a * f.a + g.a * g.a); 570 out.v = temp * (- g.a * f.v + f.a * g.v); 571 return out; 572} 573 574 575// pow -- base is a differentiatble function, exponent is a constant. 576// (a+da)^p ~= a^p + p*a^(p-1) da 577template <typename T, int N> inline 578Jet<T, N> pow(const Jet<T, N>& f, double g) { 579 Jet<T, N> out; 580 out.a = pow(f.a, g); 581 T const temp = g * pow(f.a, g - T(1.0)); 582 out.v = temp * f.v; 583 return out; 584} 585 586// pow -- base is a constant, exponent is a differentiable function. 587// (a)^(p+dp) ~= a^p + a^p log(a) dp 588template <typename T, int N> inline 589Jet<T, N> pow(double f, const Jet<T, N>& g) { 590 Jet<T, N> out; 591 out.a = pow(f, g.a); 592 T const temp = log(f) * out.a; 593 out.v = temp * g.v; 594 return out; 595} 596 597 598// pow -- both base and exponent are differentiable functions. 599// (a+da)^(b+db) ~= a^b + b * a^(b-1) da + a^b log(a) * db 600template <typename T, int N> inline 601Jet<T, N> pow(const Jet<T, N>& f, const Jet<T, N>& g) { 602 Jet<T, N> out; 603 604 T const temp1 = pow(f.a, g.a); 605 T const temp2 = g.a * pow(f.a, g.a - T(1.0)); 606 T const temp3 = temp1 * log(f.a); 607 608 out.a = temp1; 609 out.v = temp2 * f.v + temp3 * g.v; 610 return out; 611} 612 613// Define the helper functions Eigen needs to embed Jet types. 614// 615// NOTE(keir): machine_epsilon() and precision() are missing, because they don't 616// work with nested template types (e.g. where the scalar is itself templated). 617// Among other things, this means that decompositions of Jet's does not work, 618// for example 619// 620// Matrix<Jet<T, N> ... > A, x, b; 621// ... 622// A.solve(b, &x) 623// 624// does not work and will fail with a strange compiler error. 625// 626// TODO(keir): This is an Eigen 2.0 limitation that is lifted in 3.0. When we 627// switch to 3.0, also add the rest of the specialization functionality. 628template<typename T, int N> inline const Jet<T, N>& ei_conj(const Jet<T, N>& x) { return x; } // NOLINT 629template<typename T, int N> inline const Jet<T, N>& ei_real(const Jet<T, N>& x) { return x; } // NOLINT 630template<typename T, int N> inline Jet<T, N> ei_imag(const Jet<T, N>& ) { return Jet<T, N>(0.0); } // NOLINT 631template<typename T, int N> inline Jet<T, N> ei_abs (const Jet<T, N>& x) { return fabs(x); } // NOLINT 632template<typename T, int N> inline Jet<T, N> ei_abs2(const Jet<T, N>& x) { return x * x; } // NOLINT 633template<typename T, int N> inline Jet<T, N> ei_sqrt(const Jet<T, N>& x) { return sqrt(x); } // NOLINT 634template<typename T, int N> inline Jet<T, N> ei_exp (const Jet<T, N>& x) { return exp(x); } // NOLINT 635template<typename T, int N> inline Jet<T, N> ei_log (const Jet<T, N>& x) { return log(x); } // NOLINT 636template<typename T, int N> inline Jet<T, N> ei_sin (const Jet<T, N>& x) { return sin(x); } // NOLINT 637template<typename T, int N> inline Jet<T, N> ei_cos (const Jet<T, N>& x) { return cos(x); } // NOLINT 638template<typename T, int N> inline Jet<T, N> ei_pow (const Jet<T, N>& x, Jet<T, N> y) { return pow(x, y); } // NOLINT 639 640// Note: This has to be in the ceres namespace for argument dependent lookup to 641// function correctly. Otherwise statements like CHECK_LE(x, 2.0) fail with 642// strange compile errors. 643template <typename T, int N> 644inline std::ostream &operator<<(std::ostream &s, const Jet<T, N>& z) { 645 return s << "[" << z.a << " ; " << z.v.transpose() << "]"; 646} 647 648} // namespace ceres 649 650namespace Eigen { 651 652// Creating a specialization of NumTraits enables placing Jet objects inside 653// Eigen arrays, getting all the goodness of Eigen combined with autodiff. 654template<typename T, int N> 655struct NumTraits<ceres::Jet<T, N> > { 656 typedef ceres::Jet<T, N> Real; 657 typedef ceres::Jet<T, N> NonInteger; 658 typedef ceres::Jet<T, N> Nested; 659 660 static typename ceres::Jet<T, N> dummy_precision() { 661 return ceres::Jet<T, N>(1e-12); 662 } 663 664 enum { 665 IsComplex = 0, 666 IsInteger = 0, 667 IsSigned, 668 ReadCost = 1, 669 AddCost = 1, 670 // For Jet types, multiplication is more expensive than addition. 671 MulCost = 3, 672 HasFloatingPoint = 1 673 }; 674}; 675 676} // namespace Eigen 677 678#endif // CERES_PUBLIC_JET_H_ 679