jet.h revision 79397c21138f54fcff6ec067b44b847f1f7e0e98
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//   LOG(INFO) << "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  // Constructor from scalar and vector part
196  // The use of Eigen::DenseBase allows Eigen expressions
197  // to be passed in without being fully evaluated until
198  // they are assigned to v
199  template<typename Derived>
200  Jet(const T& value, const Eigen::DenseBase<Derived> &vIn)
201    : a(value),
202      v(vIn)
203  {
204  }
205
206  // Compound operators
207  Jet<T, N>& operator+=(const Jet<T, N> &y) {
208    *this = *this + y;
209    return *this;
210  }
211
212  Jet<T, N>& operator-=(const Jet<T, N> &y) {
213    *this = *this - y;
214    return *this;
215  }
216
217  Jet<T, N>& operator*=(const Jet<T, N> &y) {
218    *this = *this * y;
219    return *this;
220  }
221
222  Jet<T, N>& operator/=(const Jet<T, N> &y) {
223    *this = *this / y;
224    return *this;
225  }
226
227  // The scalar part.
228  T a;
229
230  // The infinitesimal part.
231  //
232  // Note the Eigen::DontAlign bit is needed here because this object
233  // gets allocated on the stack and as part of other arrays and
234  // structs. Forcing the right alignment there is the source of much
235  // pain and suffering. Even if that works, passing Jets around to
236  // functions by value has problems because the C++ ABI does not
237  // guarantee alignment for function arguments.
238  //
239  // Setting the DontAlign bit prevents Eigen from using SSE for the
240  // various operations on Jets. This is a small performance penalty
241  // since the AutoDiff code will still expose much of the code as
242  // statically sized loops to the compiler. But given the subtle
243  // issues that arise due to alignment, especially when dealing with
244  // multiple platforms, it seems to be a trade off worth making.
245  Eigen::Matrix<T, N, 1, Eigen::DontAlign> v;
246};
247
248// Unary +
249template<typename T, int N> inline
250Jet<T, N> const& operator+(const Jet<T, N>& f) {
251  return f;
252}
253
254// TODO(keir): Try adding __attribute__((always_inline)) to these functions to
255// see if it causes a performance increase.
256
257// Unary -
258template<typename T, int N> inline
259Jet<T, N> operator-(const Jet<T, N>&f) {
260  return Jet<T, N>(-f.a, -f.v);
261}
262
263// Binary +
264template<typename T, int N> inline
265Jet<T, N> operator+(const Jet<T, N>& f,
266                    const Jet<T, N>& g) {
267  return Jet<T, N>(f.a + g.a, f.v + g.v);
268}
269
270// Binary + with a scalar: x + s
271template<typename T, int N> inline
272Jet<T, N> operator+(const Jet<T, N>& f, T s) {
273  return Jet<T, N>(f.a + s, f.v);
274}
275
276// Binary + with a scalar: s + x
277template<typename T, int N> inline
278Jet<T, N> operator+(T s, const Jet<T, N>& f) {
279  return Jet<T, N>(f.a + s, f.v);
280}
281
282// Binary -
283template<typename T, int N> inline
284Jet<T, N> operator-(const Jet<T, N>& f,
285                    const Jet<T, N>& g) {
286  return Jet<T, N>(f.a - g.a, f.v - g.v);
287}
288
289// Binary - with a scalar: x - s
290template<typename T, int N> inline
291Jet<T, N> operator-(const Jet<T, N>& f, T s) {
292  return Jet<T, N>(f.a - s, f.v);
293}
294
295// Binary - with a scalar: s - x
296template<typename T, int N> inline
297Jet<T, N> operator-(T s, const Jet<T, N>& f) {
298  return Jet<T, N>(s - f.a, -f.v);
299}
300
301// Binary *
302template<typename T, int N> inline
303Jet<T, N> operator*(const Jet<T, N>& f,
304                    const Jet<T, N>& g) {
305  return Jet<T, N>(f.a * g.a, f.a * g.v + f.v * g.a);
306}
307
308// Binary * with a scalar: x * s
309template<typename T, int N> inline
310Jet<T, N> operator*(const Jet<T, N>& f, T s) {
311  return Jet<T, N>(f.a * s, f.v * s);
312}
313
314// Binary * with a scalar: s * x
315template<typename T, int N> inline
316Jet<T, N> operator*(T s, const Jet<T, N>& f) {
317  return Jet<T, N>(f.a * s, f.v * s);
318}
319
320// Binary /
321template<typename T, int N> inline
322Jet<T, N> operator/(const Jet<T, N>& f,
323                    const Jet<T, N>& g) {
324  // This uses:
325  //
326  //   a + u   (a + u)(b - v)   (a + u)(b - v)
327  //   ----- = -------------- = --------------
328  //   b + v   (b + v)(b - v)        b^2
329  //
330  // which holds because v*v = 0.
331  const T g_a_inverse = T(1.0) / g.a;
332  const T f_a_by_g_a = f.a * g_a_inverse;
333  return Jet<T, N>(f.a * g_a_inverse, (f.v - f_a_by_g_a * g.v) * g_a_inverse);
334}
335
336// Binary / with a scalar: s / x
337template<typename T, int N> inline
338Jet<T, N> operator/(T s, const Jet<T, N>& g) {
339  const T minus_s_g_a_inverse2 = -s / (g.a * g.a);
340  return Jet<T, N>(s / g.a, g.v * minus_s_g_a_inverse2);
341}
342
343// Binary / with a scalar: x / s
344template<typename T, int N> inline
345Jet<T, N> operator/(const Jet<T, N>& f, T s) {
346  const T s_inverse = 1.0 / s;
347  return Jet<T, N>(f.a * s_inverse, f.v * s_inverse);
348}
349
350// Binary comparison operators for both scalars and jets.
351#define CERES_DEFINE_JET_COMPARISON_OPERATOR(op) \
352template<typename T, int N> inline \
353bool operator op(const Jet<T, N>& f, const Jet<T, N>& g) { \
354  return f.a op g.a; \
355} \
356template<typename T, int N> inline \
357bool operator op(const T& s, const Jet<T, N>& g) { \
358  return s op g.a; \
359} \
360template<typename T, int N> inline \
361bool operator op(const Jet<T, N>& f, const T& s) { \
362  return f.a op s; \
363}
364CERES_DEFINE_JET_COMPARISON_OPERATOR( <  )  // NOLINT
365CERES_DEFINE_JET_COMPARISON_OPERATOR( <= )  // NOLINT
366CERES_DEFINE_JET_COMPARISON_OPERATOR( >  )  // NOLINT
367CERES_DEFINE_JET_COMPARISON_OPERATOR( >= )  // NOLINT
368CERES_DEFINE_JET_COMPARISON_OPERATOR( == )  // NOLINT
369CERES_DEFINE_JET_COMPARISON_OPERATOR( != )  // NOLINT
370#undef CERES_DEFINE_JET_COMPARISON_OPERATOR
371
372// Pull some functions from namespace std.
373//
374// This is necessary because we want to use the same name (e.g. 'sqrt') for
375// double-valued and Jet-valued functions, but we are not allowed to put
376// Jet-valued functions inside namespace std.
377//
378// TODO(keir): Switch to "using".
379inline double abs     (double x) { return std::abs(x);      }
380inline double log     (double x) { return std::log(x);      }
381inline double exp     (double x) { return std::exp(x);      }
382inline double sqrt    (double x) { return std::sqrt(x);     }
383inline double cos     (double x) { return std::cos(x);      }
384inline double acos    (double x) { return std::acos(x);     }
385inline double sin     (double x) { return std::sin(x);      }
386inline double asin    (double x) { return std::asin(x);     }
387inline double tan     (double x) { return std::tan(x);      }
388inline double atan    (double x) { return std::atan(x);     }
389inline double sinh    (double x) { return std::sinh(x);     }
390inline double cosh    (double x) { return std::cosh(x);     }
391inline double tanh    (double x) { return std::tanh(x);     }
392inline double pow  (double x, double y) { return std::pow(x, y);   }
393inline double atan2(double y, double x) { return std::atan2(y, x); }
394
395// In general, f(a + h) ~= f(a) + f'(a) h, via the chain rule.
396
397// abs(x + h) ~= x + h or -(x + h)
398template <typename T, int N> inline
399Jet<T, N> abs(const Jet<T, N>& f) {
400  return f.a < T(0.0) ? -f : f;
401}
402
403// log(a + h) ~= log(a) + h / a
404template <typename T, int N> inline
405Jet<T, N> log(const Jet<T, N>& f) {
406  const T a_inverse = T(1.0) / f.a;
407  return Jet<T, N>(log(f.a), f.v * a_inverse);
408}
409
410// exp(a + h) ~= exp(a) + exp(a) h
411template <typename T, int N> inline
412Jet<T, N> exp(const Jet<T, N>& f) {
413  const T tmp = exp(f.a);
414  return Jet<T, N>(tmp, tmp * f.v);
415}
416
417// sqrt(a + h) ~= sqrt(a) + h / (2 sqrt(a))
418template <typename T, int N> inline
419Jet<T, N> sqrt(const Jet<T, N>& f) {
420  const T tmp = sqrt(f.a);
421  const T two_a_inverse = T(1.0) / (T(2.0) * tmp);
422  return Jet<T, N>(tmp, f.v * two_a_inverse);
423}
424
425// cos(a + h) ~= cos(a) - sin(a) h
426template <typename T, int N> inline
427Jet<T, N> cos(const Jet<T, N>& f) {
428  return Jet<T, N>(cos(f.a), - sin(f.a) * f.v);
429}
430
431// acos(a + h) ~= acos(a) - 1 / sqrt(1 - a^2) h
432template <typename T, int N> inline
433Jet<T, N> acos(const Jet<T, N>& f) {
434  const T tmp = - T(1.0) / sqrt(T(1.0) - f.a * f.a);
435  return Jet<T, N>(acos(f.a), tmp * f.v);
436}
437
438// sin(a + h) ~= sin(a) + cos(a) h
439template <typename T, int N> inline
440Jet<T, N> sin(const Jet<T, N>& f) {
441  return Jet<T, N>(sin(f.a), cos(f.a) * f.v);
442}
443
444// asin(a + h) ~= asin(a) + 1 / sqrt(1 - a^2) h
445template <typename T, int N> inline
446Jet<T, N> asin(const Jet<T, N>& f) {
447  const T tmp = T(1.0) / sqrt(T(1.0) - f.a * f.a);
448  return Jet<T, N>(asin(f.a), tmp * f.v);
449}
450
451// tan(a + h) ~= tan(a) + (1 + tan(a)^2) h
452template <typename T, int N> inline
453Jet<T, N> tan(const Jet<T, N>& f) {
454  const T tan_a = tan(f.a);
455  const T tmp = T(1.0) + tan_a * tan_a;
456  return Jet<T, N>(tan_a, tmp * f.v);
457}
458
459// atan(a + h) ~= atan(a) + 1 / (1 + a^2) h
460template <typename T, int N> inline
461Jet<T, N> atan(const Jet<T, N>& f) {
462  const T tmp = T(1.0) / (T(1.0) + f.a * f.a);
463  return Jet<T, N>(atan(f.a), tmp * f.v);
464}
465
466// sinh(a + h) ~= sinh(a) + cosh(a) h
467template <typename T, int N> inline
468Jet<T, N> sinh(const Jet<T, N>& f) {
469  return Jet<T, N>(sinh(f.a), cosh(f.a) * f.v);
470}
471
472// cosh(a + h) ~= cosh(a) + sinh(a) h
473template <typename T, int N> inline
474Jet<T, N> cosh(const Jet<T, N>& f) {
475  return Jet<T, N>(cosh(f.a), sinh(f.a) * f.v);
476}
477
478// tanh(a + h) ~= tanh(a) + (1 - tanh(a)^2) h
479template <typename T, int N> inline
480Jet<T, N> tanh(const Jet<T, N>& f) {
481  const T tanh_a = tanh(f.a);
482  const T tmp = T(1.0) - tanh_a * tanh_a;
483  return Jet<T, N>(tanh_a, tmp * f.v);
484}
485
486// Jet Classification. It is not clear what the appropriate semantics are for
487// these classifications. This picks that IsFinite and isnormal are "all"
488// operations, i.e. all elements of the jet must be finite for the jet itself
489// to be finite (or normal). For IsNaN and IsInfinite, the answer is less
490// clear. This takes a "any" approach for IsNaN and IsInfinite such that if any
491// part of a jet is nan or inf, then the entire jet is nan or inf. This leads
492// to strange situations like a jet can be both IsInfinite and IsNaN, but in
493// practice the "any" semantics are the most useful for e.g. checking that
494// derivatives are sane.
495
496// The jet is finite if all parts of the jet are finite.
497template <typename T, int N> inline
498bool IsFinite(const Jet<T, N>& f) {
499  if (!IsFinite(f.a)) {
500    return false;
501  }
502  for (int i = 0; i < N; ++i) {
503    if (!IsFinite(f.v[i])) {
504      return false;
505    }
506  }
507  return true;
508}
509
510// The jet is infinite if any part of the jet is infinite.
511template <typename T, int N> inline
512bool IsInfinite(const Jet<T, N>& f) {
513  if (IsInfinite(f.a)) {
514    return true;
515  }
516  for (int i = 0; i < N; i++) {
517    if (IsInfinite(f.v[i])) {
518      return true;
519    }
520  }
521  return false;
522}
523
524// The jet is NaN if any part of the jet is NaN.
525template <typename T, int N> inline
526bool IsNaN(const Jet<T, N>& f) {
527  if (IsNaN(f.a)) {
528    return true;
529  }
530  for (int i = 0; i < N; ++i) {
531    if (IsNaN(f.v[i])) {
532      return true;
533    }
534  }
535  return false;
536}
537
538// The jet is normal if all parts of the jet are normal.
539template <typename T, int N> inline
540bool IsNormal(const Jet<T, N>& f) {
541  if (!IsNormal(f.a)) {
542    return false;
543  }
544  for (int i = 0; i < N; ++i) {
545    if (!IsNormal(f.v[i])) {
546      return false;
547    }
548  }
549  return true;
550}
551
552// atan2(b + db, a + da) ~= atan2(b, a) + (- b da + a db) / (a^2 + b^2)
553//
554// In words: the rate of change of theta is 1/r times the rate of
555// change of (x, y) in the positive angular direction.
556template <typename T, int N> inline
557Jet<T, N> atan2(const Jet<T, N>& g, const Jet<T, N>& f) {
558  // Note order of arguments:
559  //
560  //   f = a + da
561  //   g = b + db
562
563  T const tmp = T(1.0) / (f.a * f.a + g.a * g.a);
564  return Jet<T, N>(atan2(g.a, f.a), tmp * (- g.a * f.v + f.a * g.v));
565}
566
567
568// pow -- base is a differentiable function, exponent is a constant.
569// (a+da)^p ~= a^p + p*a^(p-1) da
570template <typename T, int N> inline
571Jet<T, N> pow(const Jet<T, N>& f, double g) {
572  T const tmp = g * pow(f.a, g - T(1.0));
573  return Jet<T, N>(pow(f.a, g), tmp * f.v);
574}
575
576// pow -- base is a constant, exponent is a differentiable function.
577// (a)^(p+dp) ~= a^p + a^p log(a) dp
578template <typename T, int N> inline
579Jet<T, N> pow(double f, const Jet<T, N>& g) {
580  T const tmp = pow(f, g.a);
581  return Jet<T, N>(tmp, log(f) * tmp * g.v);
582}
583
584
585// pow -- both base and exponent are differentiable functions.
586// (a+da)^(b+db) ~= a^b + b * a^(b-1) da + a^b log(a) * db
587template <typename T, int N> inline
588Jet<T, N> pow(const Jet<T, N>& f, const Jet<T, N>& g) {
589  T const tmp1 = pow(f.a, g.a);
590  T const tmp2 = g.a * pow(f.a, g.a - T(1.0));
591  T const tmp3 = tmp1 * log(f.a);
592
593  return Jet<T, N>(tmp1, tmp2 * f.v + tmp3 * g.v);
594}
595
596// Define the helper functions Eigen needs to embed Jet types.
597//
598// NOTE(keir): machine_epsilon() and precision() are missing, because they don't
599// work with nested template types (e.g. where the scalar is itself templated).
600// Among other things, this means that decompositions of Jet's does not work,
601// for example
602//
603//   Matrix<Jet<T, N> ... > A, x, b;
604//   ...
605//   A.solve(b, &x)
606//
607// does not work and will fail with a strange compiler error.
608//
609// TODO(keir): This is an Eigen 2.0 limitation that is lifted in 3.0. When we
610// switch to 3.0, also add the rest of the specialization functionality.
611template<typename T, int N> inline const Jet<T, N>& ei_conj(const Jet<T, N>& x) { return x;              }  // NOLINT
612template<typename T, int N> inline const Jet<T, N>& ei_real(const Jet<T, N>& x) { return x;              }  // NOLINT
613template<typename T, int N> inline       Jet<T, N>  ei_imag(const Jet<T, N>&  ) { return Jet<T, N>(0.0); }  // NOLINT
614template<typename T, int N> inline       Jet<T, N>  ei_abs (const Jet<T, N>& x) { return fabs(x);        }  // NOLINT
615template<typename T, int N> inline       Jet<T, N>  ei_abs2(const Jet<T, N>& x) { return x * x;          }  // NOLINT
616template<typename T, int N> inline       Jet<T, N>  ei_sqrt(const Jet<T, N>& x) { return sqrt(x);        }  // NOLINT
617template<typename T, int N> inline       Jet<T, N>  ei_exp (const Jet<T, N>& x) { return exp(x);         }  // NOLINT
618template<typename T, int N> inline       Jet<T, N>  ei_log (const Jet<T, N>& x) { return log(x);         }  // NOLINT
619template<typename T, int N> inline       Jet<T, N>  ei_sin (const Jet<T, N>& x) { return sin(x);         }  // NOLINT
620template<typename T, int N> inline       Jet<T, N>  ei_cos (const Jet<T, N>& x) { return cos(x);         }  // NOLINT
621template<typename T, int N> inline       Jet<T, N>  ei_tan (const Jet<T, N>& x) { return tan(x);         }  // NOLINT
622template<typename T, int N> inline       Jet<T, N>  ei_atan(const Jet<T, N>& x) { return atan(x);        }  // NOLINT
623template<typename T, int N> inline       Jet<T, N>  ei_sinh(const Jet<T, N>& x) { return sinh(x);        }  // NOLINT
624template<typename T, int N> inline       Jet<T, N>  ei_cosh(const Jet<T, N>& x) { return cosh(x);        }  // NOLINT
625template<typename T, int N> inline       Jet<T, N>  ei_tanh(const Jet<T, N>& x) { return tanh(x);        }  // NOLINT
626template<typename T, int N> inline       Jet<T, N>  ei_pow (const Jet<T, N>& x, Jet<T, N> y) { return pow(x, y); }  // NOLINT
627
628// Note: This has to be in the ceres namespace for argument dependent lookup to
629// function correctly. Otherwise statements like CHECK_LE(x, 2.0) fail with
630// strange compile errors.
631template <typename T, int N>
632inline std::ostream &operator<<(std::ostream &s, const Jet<T, N>& z) {
633  return s << "[" << z.a << " ; " << z.v.transpose() << "]";
634}
635
636}  // namespace ceres
637
638namespace Eigen {
639
640// Creating a specialization of NumTraits enables placing Jet objects inside
641// Eigen arrays, getting all the goodness of Eigen combined with autodiff.
642template<typename T, int N>
643struct NumTraits<ceres::Jet<T, N> > {
644  typedef ceres::Jet<T, N> Real;
645  typedef ceres::Jet<T, N> NonInteger;
646  typedef ceres::Jet<T, N> Nested;
647
648  static typename ceres::Jet<T, N> dummy_precision() {
649    return ceres::Jet<T, N>(1e-12);
650  }
651
652  static inline Real epsilon() { return Real(std::numeric_limits<T>::epsilon()); }
653
654  enum {
655    IsComplex = 0,
656    IsInteger = 0,
657    IsSigned,
658    ReadCost = 1,
659    AddCost = 1,
660    // For Jet types, multiplication is more expensive than addition.
661    MulCost = 3,
662    HasFloatingPoint = 1,
663    RequireInitialization = 1
664  };
665};
666
667}  // namespace Eigen
668
669#endif  // CERES_PUBLIC_JET_H_
670