NumericalDiff.h
1 // -*- coding: utf-8
2 // vim: set fileencoding=utf-8
3 
4 // This file is part of Eigen, a lightweight C++ template library
5 // for linear algebra.
6 //
7 // Copyright (C) 2009 Thomas Capricelli <orzel@freehackers.org>
8 //
9 // This Source Code Form is subject to the terms of the Mozilla
10 // Public License v. 2.0. If a copy of the MPL was not distributed
11 // with this file, You can obtain one at http://mozilla.org/MPL/2.0/.
12 
13 #ifndef EIGEN_NUMERICAL_DIFF_H
14 #define EIGEN_NUMERICAL_DIFF_H
15 
16 namespace Eigen {
17 
18 enum NumericalDiffMode {
19  Forward,
20  Central
21 };
22 
23 
35 template<typename _Functor, NumericalDiffMode mode=Forward>
36 class NumericalDiff : public _Functor
37 {
38 public:
39  typedef _Functor Functor;
40  typedef typename Functor::Scalar Scalar;
41  typedef typename Functor::InputType InputType;
42  typedef typename Functor::ValueType ValueType;
43  typedef typename Functor::JacobianType JacobianType;
44 
45  NumericalDiff(Scalar _epsfcn=0.) : Functor(), epsfcn(_epsfcn) {}
46  NumericalDiff(const Functor& f, Scalar _epsfcn=0.) : Functor(f), epsfcn(_epsfcn) {}
47 
48  // forward constructors
49  template<typename T0>
50  NumericalDiff(const T0& a0) : Functor(a0), epsfcn(0) {}
51  template<typename T0, typename T1>
52  NumericalDiff(const T0& a0, const T1& a1) : Functor(a0, a1), epsfcn(0) {}
53  template<typename T0, typename T1, typename T2>
54  NumericalDiff(const T0& a0, const T1& a1, const T2& a2) : Functor(a0, a1, a2), epsfcn(0) {}
55 
56  enum {
57  InputsAtCompileTime = Functor::InputsAtCompileTime,
58  ValuesAtCompileTime = Functor::ValuesAtCompileTime
59  };
60 
64  int df(const InputType& _x, JacobianType &jac) const
65  {
66  /* Local variables */
67  Scalar h;
68  int nfev=0;
69  const typename InputType::Index n = _x.size();
70  const Scalar eps = internal::sqrt(((std::max)(epsfcn,NumTraits<Scalar>::epsilon() )));
71  ValueType val1, val2;
72  InputType x = _x;
73  // TODO : we should do this only if the size is not already known
74  val1.resize(Functor::values());
75  val2.resize(Functor::values());
76 
77  // initialization
78  switch(mode) {
79  case Forward:
80  // compute f(x)
81  Functor::operator()(x, val1); nfev++;
82  break;
83  case Central:
84  // do nothing
85  break;
86  default:
87  assert(false);
88  };
89 
90  // Function Body
91  for (int j = 0; j < n; ++j) {
92  h = eps * internal::abs(x[j]);
93  if (h == 0.) {
94  h = eps;
95  }
96  switch(mode) {
97  case Forward:
98  x[j] += h;
99  Functor::operator()(x, val2);
100  nfev++;
101  x[j] = _x[j];
102  jac.col(j) = (val2-val1)/h;
103  break;
104  case Central:
105  x[j] += h;
106  Functor::operator()(x, val2); nfev++;
107  x[j] -= 2*h;
108  Functor::operator()(x, val1); nfev++;
109  x[j] = _x[j];
110  jac.col(j) = (val2-val1)/(2*h);
111  break;
112  default:
113  assert(false);
114  };
115  }
116  return nfev;
117  }
118 private:
119  Scalar epsfcn;
120 
121  NumericalDiff& operator=(const NumericalDiff&);
122 };
123 
124 } // end namespace Eigen
125 
126 //vim: ai ts=4 sts=4 et sw=4
127 #endif // EIGEN_NUMERICAL_DIFF_H
128