Fix compile errors on VC++ 2012 RC1.
[blender.git] / extern / libmv / third_party / ceres / include / ceres / jet.h
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
166 // Visual Studio 2012 or older version
167 #if defined(_MSC_VER) && _MSC_VER <= 1700
168 namespace std {
169 inline bool isfinite(double x) { return _finite(x);                }
170 inline bool isinf   (double x) { return !_finite(x) && !_isnan(x); }
171 inline bool isnan   (double x) { return _isnan(x);                 }
172 inline bool isnormal(double x) { return _finite(x) && x != 0.0;    }
173 }  // namespace std
174 #endif
175
176 namespace ceres {
177
178 template <typename T, int N>
179 struct Jet {
180   enum { DIMENSION = N };
181
182   // Default-construct "a" because otherwise this can lead to false errors about
183   // uninitialized uses when other classes relying on default constructed T
184   // (where T is a Jet<T, N>). This usually only happens in opt mode. Note that
185   // the C++ standard mandates that e.g. default constructed doubles are
186   // initialized to 0.0; see sections 8.5 of the C++03 standard.
187   Jet() : a() {}
188
189   // Constructor from scalar: a + 0.
190   explicit Jet(const T& value) {
191     a = value;
192     v.setZero();
193   }
194
195   // Constructor from scalar plus variable: a + t_i.
196   Jet(const T& value, int k) {
197     a = value;
198     v.setZero();
199     v[k] = T(1.0);
200   }
201
202   /*
203
204   // Construct from an array where the first element is the scalar.
205   // This is templated to support converting from other data types.
206   template<typename D>
207   Jet(const D* scalar_and_derivatives) {
208     a = T(scalar_and_derivatives[0]);
209     v = Eigen::Map<const Eigen::Matrix<D, N, 1> >(
210         scalar_and_derivatives + 1, N).cast<T>();
211   }
212   */
213
214   // Compound operators
215   Jet<T, N>& operator+=(const Jet<T, N> &y) {
216     *this = *this + y;
217     return *this;
218   }
219
220   Jet<T, N>& operator-=(const Jet<T, N> &y) {
221     *this = *this - y;
222     return *this;
223   }
224
225   Jet<T, N>& operator*=(const Jet<T, N> &y) {
226     *this = *this * y;
227     return *this;
228   }
229
230   Jet<T, N>& operator/=(const Jet<T, N> &y) {
231     *this = *this / y;
232     return *this;
233   }
234
235   T a;  // The scalar part.
236   Eigen::Matrix<T, N, 1> v;  // The infinitesimal part.
237 };
238
239 // Unary +
240 template<typename T, int N> inline
241 Jet<T, N> const& operator+(const Jet<T, N>& f) {
242   return f;
243 }
244
245 // TODO(keir): Try adding __attribute__((always_inline)) to these functions to
246 // see if it causes a performance increase.
247
248 // Unary -
249 template<typename T, int N> inline
250 Jet<T, N> operator-(const Jet<T, N>&f) {
251   Jet<T, N> g;
252   g.a = -f.a;
253   g.v = -f.v;
254   return g;
255 }
256
257 // Binary +
258 template<typename T, int N> inline
259 Jet<T, N> operator+(const Jet<T, N>& f,
260                     const Jet<T, N>& g) {
261   Jet<T, N> h;
262   h.a = f.a + g.a;
263   h.v = f.v + g.v;
264   return h;
265 }
266
267 // Binary + with a scalar: x + s
268 template<typename T, int N> inline
269 Jet<T, N> operator+(const Jet<T, N>& f, T s) {
270   Jet<T, N> h;
271   h.a = f.a + s;
272   h.v = f.v;
273   return h;
274 }
275
276 // Binary + with a scalar: s + x
277 template<typename T, int N> inline
278 Jet<T, N> operator+(T s, const Jet<T, N>& f) {
279   Jet<T, N> h;
280   h.a = f.a + s;
281   h.v = f.v;
282   return h;
283 }
284
285 // Binary -
286 template<typename T, int N> inline
287 Jet<T, N> operator-(const Jet<T, N>& f,
288                     const Jet<T, N>& g) {
289   Jet<T, N> h;
290   h.a = f.a - g.a;
291   h.v = f.v - g.v;
292   return h;
293 }
294
295 // Binary - with a scalar: x - s
296 template<typename T, int N> inline
297 Jet<T, N> operator-(const Jet<T, N>& f, T s) {
298   Jet<T, N> h;
299   h.a = f.a - s;
300   h.v = f.v;
301   return h;
302 }
303
304 // Binary - with a scalar: s - x
305 template<typename T, int N> inline
306 Jet<T, N> operator-(T s, const Jet<T, N>& f) {
307   Jet<T, N> h;
308   h.a = s - f.a;
309   h.v = -f.v;
310   return h;
311 }
312
313 // Binary *
314 template<typename T, int N> inline
315 Jet<T, N> operator*(const Jet<T, N>& f,
316                     const Jet<T, N>& g) {
317   Jet<T, N> h;
318   h.a = f.a * g.a;
319   h.v = f.a * g.v + f.v * g.a;
320   return h;
321 }
322
323 // Binary * with a scalar: x * s
324 template<typename T, int N> inline
325 Jet<T, N> operator*(const Jet<T, N>& f, T s) {
326   Jet<T, N> h;
327   h.a = f.a * s;
328   h.v = f.v * s;
329   return h;
330 }
331
332 // Binary * with a scalar: s * x
333 template<typename T, int N> inline
334 Jet<T, N> operator*(T s, const Jet<T, N>& f) {
335   Jet<T, N> h;
336   h.a = f.a * s;
337   h.v = f.v * s;
338   return h;
339 }
340
341 // Binary /
342 template<typename T, int N> inline
343 Jet<T, N> operator/(const Jet<T, N>& f,
344                     const Jet<T, N>& g) {
345   Jet<T, N> h;
346   // This uses:
347   //
348   //   a + u   (a + u)(b - v)   (a + u)(b - v)
349   //   ----- = -------------- = --------------
350   //   b + v   (b + v)(b - v)        b^2
351   //
352   // which holds because v*v = 0.
353   h.a = f.a / g.a;
354   h.v = (f.v - f.a / g.a * g.v) / g.a;
355   return h;
356 }
357
358 // Binary / with a scalar: s / x
359 template<typename T, int N> inline
360 Jet<T, N> operator/(T s, const Jet<T, N>& g) {
361   Jet<T, N> h;
362   h.a = s / g.a;
363   h.v = - s * g.v / (g.a * g.a);
364   return h;
365 }
366
367 // Binary / with a scalar: x / s
368 template<typename T, int N> inline
369 Jet<T, N> operator/(const Jet<T, N>& f, T s) {
370   Jet<T, N> h;
371   h.a = f.a / s;
372   h.v = f.v / s;
373   return h;
374 }
375
376 // Binary comparison operators for both scalars and jets.
377 #define CERES_DEFINE_JET_COMPARISON_OPERATOR(op) \
378 template<typename T, int N> inline \
379 bool operator op(const Jet<T, N>& f, const Jet<T, N>& g) { \
380   return f.a op g.a; \
381 } \
382 template<typename T, int N> inline \
383 bool operator op(const T& s, const Jet<T, N>& g) { \
384   return s op g.a; \
385 } \
386 template<typename T, int N> inline \
387 bool operator op(const Jet<T, N>& f, const T& s) { \
388   return f.a op s; \
389 }
390 CERES_DEFINE_JET_COMPARISON_OPERATOR( <  )  // NOLINT
391 CERES_DEFINE_JET_COMPARISON_OPERATOR( <= )  // NOLINT
392 CERES_DEFINE_JET_COMPARISON_OPERATOR( >  )  // NOLINT
393 CERES_DEFINE_JET_COMPARISON_OPERATOR( >= )  // NOLINT
394 CERES_DEFINE_JET_COMPARISON_OPERATOR( == )  // NOLINT
395 CERES_DEFINE_JET_COMPARISON_OPERATOR( != )  // NOLINT
396 #undef CERES_DEFINE_JET_COMPARISON_OPERATOR
397
398 // Pull some functions from namespace std.
399 //
400 // This is necessary because we want to use the same name (e.g. 'sqrt') for
401 // double-valued and Jet-valued functions, but we are not allowed to put
402 // Jet-valued functions inside namespace std.
403 //
404 // Missing: cosh, sinh, tanh, tan
405 // TODO(keir): Switch to "using".
406 inline double abs     (double x) { return std::abs(x);      }
407 inline double log     (double x) { return std::log(x);      }
408 inline double exp     (double x) { return std::exp(x);      }
409 inline double sqrt    (double x) { return std::sqrt(x);     }
410 inline double cos     (double x) { return std::cos(x);      }
411 inline double acos    (double x) { return std::acos(x);     }
412 inline double sin     (double x) { return std::sin(x);      }
413 inline double asin    (double x) { return std::asin(x);     }
414 inline bool   isfinite(double x) { return std::isfinite(x); }
415 inline bool   isinf   (double x) { return std::isinf(x);    }
416 inline bool   isnan   (double x) { return std::isnan(x);    }
417 inline bool   isnormal(double x) { return std::isnormal(x); }
418 inline double pow  (double x, double y) { return std::pow(x, y);   }
419 inline double atan2(double y, double x) { return std::atan2(y, x); }
420
421 // In general, f(a + h) ~= f(a) + f'(a) h, via the chain rule.
422
423 // abs(x + h) ~= x + h or -(x + h)
424 template <typename T, int N> inline
425 Jet<T, N> abs(const Jet<T, N>& f) {
426   return f.a < T(0.0) ? -f : f;
427 }
428
429 // log(a + h) ~= log(a) + h / a
430 template <typename T, int N> inline
431 Jet<T, N> log(const Jet<T, N>& f) {
432   Jet<T, N> g;
433   g.a = log(f.a);
434   g.v = f.v / f.a;
435   return g;
436 }
437
438 // exp(a + h) ~= exp(a) + exp(a) h
439 template <typename T, int N> inline
440 Jet<T, N> exp(const Jet<T, N>& f) {
441   Jet<T, N> g;
442   g.a = exp(f.a);
443   g.v = g.a * f.v;
444   return g;
445 }
446
447 // sqrt(a + h) ~= sqrt(a) + h / (2 sqrt(a))
448 template <typename T, int N> inline
449 Jet<T, N> sqrt(const Jet<T, N>& f) {
450   Jet<T, N> g;
451   g.a = sqrt(f.a);
452   g.v = f.v / (T(2.0) * g.a);
453   return g;
454 }
455
456 // cos(a + h) ~= cos(a) - sin(a) h
457 template <typename T, int N> inline
458 Jet<T, N> cos(const Jet<T, N>& f) {
459   Jet<T, N> g;
460   g.a = cos(f.a);
461   T sin_a = sin(f.a);
462   g.v = - sin_a * f.v;
463   return g;
464 }
465
466 // acos(a + h) ~= acos(a) - 1 / sqrt(1 - a^2) h
467 template <typename T, int N> inline
468 Jet<T, N> acos(const Jet<T, N>& f) {
469   Jet<T, N> g;
470   g.a = acos(f.a);
471   g.v = - T(1.0) / sqrt(T(1.0) - f.a * f.a) * f.v;
472   return g;
473 }
474
475 // sin(a + h) ~= sin(a) + cos(a) h
476 template <typename T, int N> inline
477 Jet<T, N> sin(const Jet<T, N>& f) {
478   Jet<T, N> g;
479   g.a = sin(f.a);
480   T cos_a = cos(f.a);
481   g.v = cos_a * f.v;
482   return g;
483 }
484
485 // asin(a + h) ~= asin(a) + 1 / sqrt(1 - a^2) h
486 template <typename T, int N> inline
487 Jet<T, N> asin(const Jet<T, N>& f) {
488   Jet<T, N> g;
489   g.a = asin(f.a);
490   g.v = T(1.0) / sqrt(T(1.0) - f.a * f.a) * f.v;
491   return g;
492 }
493
494 // Jet Classification. It is not clear what the appropriate semantics are for
495 // these classifications. This picks that isfinite and isnormal are "all"
496 // operations, i.e. all elements of the jet must be finite for the jet itself to
497 // be finite (or normal). For isnan and isinf, the answer is less clear. This
498 // takes a "any" approach for isnan and isinf such that if any part of a jet is
499 // nan or inf, then the entire jet is nan or inf. This leads to strange
500 // situations like a jet can be both isinf and isnan, but in practice the "any"
501 // semantics are the most useful for e.g. checking that derivatives are sane.
502
503 // The jet is finite if all parts of the jet are finite.
504 template <typename T, int N> inline
505 bool isfinite(const Jet<T, N>& f) {
506   if (!isfinite(f.a)) {
507     return false;
508   }
509   for (int i = 0; i < N; ++i) {
510     if (!isfinite(f.v[i])) {
511       return false;
512     }
513   }
514   return true;
515 }
516
517 // The jet is infinite if any part of the jet is infinite.
518 template <typename T, int N> inline
519 bool isinf(const Jet<T, N>& f) {
520   if (isinf(f.a)) {
521     return true;
522   }
523   for (int i = 0; i < N; i++) {
524     if (isinf(f.v[i])) {
525       return true;
526     }
527   }
528   return false;
529 }
530
531 // The jet is NaN if any part of the jet is NaN.
532 template <typename T, int N> inline
533 bool isnan(const Jet<T, N>& f) {
534   if (isnan(f.a)) {
535     return true;
536   }
537   for (int i = 0; i < N; ++i) {
538     if (isnan(f.v[i])) {
539       return true;
540     }
541   }
542   return false;
543 }
544
545 // The jet is normal if all parts of the jet are normal.
546 template <typename T, int N> inline
547 bool isnormal(const Jet<T, N>& f) {
548   if (!isnormal(f.a)) {
549     return false;
550   }
551   for (int i = 0; i < N; ++i) {
552     if (!isnormal(f.v[i])) {
553       return false;
554     }
555   }
556   return true;
557 }
558
559 // atan2(b + db, a + da) ~= atan2(b, a) + (- b da + a db) / (a^2 + b^2)
560 //
561 // In words: the rate of change of theta is 1/r times the rate of
562 // change of (x, y) in the positive angular direction.
563 template <typename T, int N> inline
564 Jet<T, N> atan2(const Jet<T, N>& g, const Jet<T, N>& f) {
565   // Note order of arguments:
566   //
567   //   f = a + da
568   //   g = b + db
569
570   Jet<T, N> out;
571
572   out.a = atan2(g.a, f.a);
573
574   T const temp = T(1.0) / (f.a * f.a + g.a * g.a);
575   out.v = temp * (- g.a * f.v + f.a * g.v);
576   return out;
577 }
578
579
580 // pow -- base is a differentiatble function, exponent is a constant.
581 // (a+da)^p ~= a^p + p*a^(p-1) da
582 template <typename T, int N> inline
583 Jet<T, N> pow(const Jet<T, N>& f, double g) {
584   Jet<T, N> out;
585   out.a = pow(f.a, g);
586   T const temp = g * pow(f.a, g - T(1.0));
587   out.v = temp * f.v;
588   return out;
589 }
590
591 // pow -- base is a constant, exponent is a differentiable function.
592 // (a)^(p+dp) ~= a^p + a^p log(a) dp
593 template <typename T, int N> inline
594 Jet<T, N> pow(double f, const Jet<T, N>& g) {
595   Jet<T, N> out;
596   out.a = pow(f, g.a);
597   T const temp = log(f) * out.a;
598   out.v = temp * g.v;
599   return out;
600 }
601
602
603 // pow -- both base and exponent are differentiable functions.
604 // (a+da)^(b+db) ~= a^b + b * a^(b-1) da + a^b log(a) * db
605 template <typename T, int N> inline
606 Jet<T, N> pow(const Jet<T, N>& f, const Jet<T, N>& g) {
607   Jet<T, N> out;
608
609   T const temp1 = pow(f.a, g.a);
610   T const temp2 = g.a * pow(f.a, g.a - T(1.0));
611   T const temp3 = temp1 * log(f.a);
612
613   out.a = temp1;
614   out.v = temp2 * f.v + temp3 * g.v;
615   return out;
616 }
617
618 // Define the helper functions Eigen needs to embed Jet types.
619 //
620 // NOTE(keir): machine_epsilon() and precision() are missing, because they don't
621 // work with nested template types (e.g. where the scalar is itself templated).
622 // Among other things, this means that decompositions of Jet's does not work,
623 // for example
624 //
625 //   Matrix<Jet<T, N> ... > A, x, b;
626 //   ...
627 //   A.solve(b, &x)
628 //
629 // does not work and will fail with a strange compiler error.
630 //
631 // TODO(keir): This is an Eigen 2.0 limitation that is lifted in 3.0. When we
632 // switch to 3.0, also add the rest of the specialization functionality.
633 template<typename T, int N> inline const Jet<T, N>& ei_conj(const Jet<T, N>& x) { return x;              }  // NOLINT
634 template<typename T, int N> inline const Jet<T, N>& ei_real(const Jet<T, N>& x) { return x;              }  // NOLINT
635 template<typename T, int N> inline       Jet<T, N>  ei_imag(const Jet<T, N>&  ) { return Jet<T, N>(0.0); }  // NOLINT
636 template<typename T, int N> inline       Jet<T, N>  ei_abs (const Jet<T, N>& x) { return fabs(x);        }  // NOLINT
637 template<typename T, int N> inline       Jet<T, N>  ei_abs2(const Jet<T, N>& x) { return x * x;          }  // NOLINT
638 template<typename T, int N> inline       Jet<T, N>  ei_sqrt(const Jet<T, N>& x) { return sqrt(x);        }  // NOLINT
639 template<typename T, int N> inline       Jet<T, N>  ei_exp (const Jet<T, N>& x) { return exp(x);         }  // NOLINT
640 template<typename T, int N> inline       Jet<T, N>  ei_log (const Jet<T, N>& x) { return log(x);         }  // NOLINT
641 template<typename T, int N> inline       Jet<T, N>  ei_sin (const Jet<T, N>& x) { return sin(x);         }  // NOLINT
642 template<typename T, int N> inline       Jet<T, N>  ei_cos (const Jet<T, N>& x) { return cos(x);         }  // NOLINT
643 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
644
645 // Note: This has to be in the ceres namespace for argument dependent lookup to
646 // function correctly. Otherwise statements like CHECK_LE(x, 2.0) fail with
647 // strange compile errors.
648 template <typename T, int N>
649 inline std::ostream &operator<<(std::ostream &s, const Jet<T, N>& z) {
650   return s << "[" << z.a << " ; " << z.v.transpose() << "]";
651 }
652
653 // A jet traits class to make it easier to work with mixed auto / numeric diff.
654 template<typename T>
655 struct JetOps {
656   static bool IsScalar() {
657     return true;
658   }
659   static T GetScalar(const T& t) {
660     return t;
661   }
662   static void SetScalar(const T& scalar, T* t) {
663     *t = scalar;
664   }
665   static void ScaleDerivative(double scale_by, T *value) {
666     // For double, there is no derivative to scale.
667   }
668 };
669
670 template<typename T, int N>
671 struct JetOps<Jet<T, N> > {
672   static bool IsScalar() {
673     return false;
674   }
675   static T GetScalar(const Jet<T, N>& t) {
676     return t.a;
677   }
678   static void SetScalar(const T& scalar, Jet<T, N>* t) {
679     t->a = scalar;
680   }
681   static void ScaleDerivative(double scale_by, Jet<T, N> *value) {
682     value->v *= scale_by;
683   }
684 };
685
686 template<typename FunctionType, int kNumArgs, typename ArgumentType>
687 struct Chain {
688   static ArgumentType Rule(const FunctionType &f,
689                            const FunctionType dfdx[kNumArgs],
690                            const ArgumentType x[kNumArgs]) {
691     // In the default case of scalars, there's nothing to do since there are no
692     // derivatives to propagate. 
693     return f;
694   }
695 };
696
697 // XXX Add documentation here!
698 template<typename FunctionType, int kNumArgs, typename T, int N>
699 struct Chain<FunctionType, kNumArgs, Jet<T, N> > {
700   static Jet<T, N> Rule(const FunctionType &f,
701                         const FunctionType dfdx[kNumArgs],
702                         const Jet<T, N> x[kNumArgs]) {
703     // x is itself a function of another variable ("z"); what this function
704     // needs to return is "f", but with the derivative with respect to z
705     // attached to the jet. So combine the derivative part of x's jets to form
706     // a Jacobian matrix between x and z (i.e. dx/dz).
707     Eigen::Matrix<T, kNumArgs, N> dxdz;
708     for (int i = 0; i < kNumArgs; ++i) {
709       dxdz.row(i) = x[i].v.transpose();
710     }
711
712     // Map the input gradient dfdx into an Eigen row vector.
713     Eigen::Map<const Eigen::Matrix<FunctionType, 1, kNumArgs> >
714         vector_dfdx(dfdx, 1, kNumArgs);
715
716     // Now apply the chain rule to obtain df/dz. Combine the derivative with
717     // the scalar part to obtain f with full derivative information.
718     Jet<T, N> jet_f;
719     jet_f.a = f;
720     jet_f.v = vector_dfdx.template cast<T>() * dxdz;  // Also known as dfdz.
721     return jet_f;
722   }
723 };
724
725 }  // namespace ceres
726
727 namespace Eigen {
728
729 // Creating a specialization of NumTraits enables placing Jet objects inside
730 // Eigen arrays, getting all the goodness of Eigen combined with autodiff.
731 template<typename T, int N>
732 struct NumTraits<ceres::Jet<T, N> > {
733   typedef ceres::Jet<T, N> Real;
734   typedef ceres::Jet<T, N> NonInteger;
735   typedef ceres::Jet<T, N> Nested;
736
737   static typename ceres::Jet<T, N> dummy_precision() {
738     return ceres::Jet<T, N>(1e-12);
739   }
740
741   enum {
742     IsComplex = 0,
743     IsInteger = 0,
744     IsSigned,
745     ReadCost = 1,
746     AddCost = 1,
747     // For Jet types, multiplication is more expensive than addition.
748     MulCost = 3,
749     HasFloatingPoint = 1
750   };
751 };
752
753 }  // namespace Eigen
754
755 #endif  // CERES_PUBLIC_JET_H_