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


gtsam
Author(s):
autogenerated on Sat May 8 2021 02:42:26