Blender  V3.3
function_derivative.h
Go to the documentation of this file.
1 // Copyright (c) 2007, 2008, 2009 libmv authors.
2 //
3 // Permission is hereby granted, free of charge, to any person obtaining a copy
4 // of this software and associated documentation files (the "Software"), to
5 // deal in the Software without restriction, including without limitation the
6 // rights to use, copy, modify, merge, publish, distribute, sublicense, and/or
7 // sell copies of the Software, and to permit persons to whom the Software is
8 // furnished to do so, subject to the following conditions:
9 //
10 // The above copyright notice and this permission notice shall be included in
11 // all copies or substantial portions of the Software.
12 //
13 // THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, EXPRESS OR
14 // IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF MERCHANTABILITY,
15 // FITNESS FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT. IN NO EVENT SHALL THE
16 // AUTHORS OR COPYRIGHT HOLDERS BE LIABLE FOR ANY CLAIM, DAMAGES OR OTHER
17 // LIABILITY, WHETHER IN AN ACTION OF CONTRACT, TORT OR OTHERWISE, ARISING
18 // FROM, OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR OTHER DEALINGS
19 // IN THE SOFTWARE.
20 
21 #ifndef LIBMV_NUMERIC_DERIVATIVE_H
22 #define LIBMV_NUMERIC_DERIVATIVE_H
23 
24 #include <cmath>
25 
26 #include "libmv/logging/logging.h"
27 #include "libmv/numeric/numeric.h"
28 
29 namespace libmv {
30 
31 // Numeric derivative of a function.
32 // TODO(keir): Consider adding a quadratic approximation.
33 
37 };
38 
39 template <typename Function, NumericJacobianMode mode = CENTRAL>
41  public:
42  typedef typename Function::XMatrixType Parameters;
43  typedef typename Function::XMatrixType::RealScalar XScalar;
44  typedef typename Function::FMatrixType FMatrixType;
45  typedef Matrix<typename Function::FMatrixType::RealScalar,
46  Function::FMatrixType::RowsAtCompileTime,
47  Function::XMatrixType::RowsAtCompileTime>
49 
50  NumericJacobian(const Function& f) : f_(f) {}
51 
52  // TODO(keir): Perhaps passing the jacobian back by value is not a good idea.
54  // Empirically determined constant.
55  Parameters eps = x.array().abs() * XScalar(1e-5);
56  // To handle cases where a paremeter is exactly zero, instead use the mean
57  // eps for the other dimensions.
58  XScalar mean_eps = eps.sum() / eps.rows();
59  if (mean_eps == XScalar(0)) {
60  // TODO(keir): Do something better here.
61  mean_eps = 1e-8; // ~sqrt(machine precision).
62  }
63  // TODO(keir): Elimininate this needless function evaluation for the
64  // central difference case.
65  FMatrixType fx = f_(x);
66  const int rows = fx.rows();
67  const int cols = x.rows();
68  JMatrixType jacobian(rows, cols);
69  Parameters x_plus_delta = x;
70  for (int c = 0; c < cols; ++c) {
71  if (eps(c) == XScalar(0)) {
72  eps(c) = mean_eps;
73  }
74  x_plus_delta(c) = x(c) + eps(c);
75  jacobian.col(c) = f_(x_plus_delta);
76 
77  XScalar one_over_h = 1 / eps(c);
78  if (mode == CENTRAL) {
79  x_plus_delta(c) = x(c) - eps(c);
80  jacobian.col(c) -= f_(x_plus_delta);
81  one_over_h /= 2;
82  } else {
83  jacobian.col(c) -= fx;
84  }
85  x_plus_delta(c) = x(c);
86  jacobian.col(c) = jacobian.col(c) * one_over_h;
87  }
88  return jacobian;
89  }
90 
91  private:
92  const Function& f_;
93 };
94 
95 template <typename Function, typename Jacobian>
96 bool CheckJacobian(const Function& f, const typename Function::XMatrixType& x) {
97  Jacobian j_analytic(f);
98  NumericJacobian<Function> j_numeric(f);
99 
100  typename NumericJacobian<Function>::JMatrixType J_numeric = j_numeric(x);
101  typename NumericJacobian<Function>::JMatrixType J_analytic = j_analytic(x);
102  LG << J_numeric - J_analytic;
103  return true;
104 }
105 
106 } // namespace libmv
107 
108 #endif // LIBMV_NUMERIC_DERIVATIVE_H
ATTR_WARN_UNUSED_RESULT const BMVert const BMEdge * e
JMatrixType operator()(const Parameters &x)
Function::XMatrixType Parameters
Function::XMatrixType::RealScalar XScalar
Function::FMatrixType FMatrixType
NumericJacobian(const Function &f)
Matrix< typename Function::FMatrixType::RealScalar, Function::FMatrixType::RowsAtCompileTime, Function::XMatrixType::RowsAtCompileTime > JMatrixType
#define LG
static unsigned c
Definition: RandGen.cpp:83
bool CheckJacobian(const Function &f, const typename Function::XMatrixType &x)
const btScalar eps
Definition: poly34.cpp:11