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 // Eigen is free software; you can redistribute it and/or
10 // modify it under the terms of the GNU Lesser General Public
11 // License as published by the Free Software Foundation; either
12 // version 3 of the License, or (at your option) any later version.
13 //
14 // Alternatively, you can redistribute it and/or
15 // modify it under the terms of the GNU General Public License as
16 // published by the Free Software Foundation; either version 2 of
17 // the License, or (at your option) any later version.
18 //
19 // Eigen is distributed in the hope that it will be useful, but WITHOUT ANY
20 // WARRANTY; without even the implied warranty of MERCHANTABILITY or FITNESS
21 // FOR A PARTICULAR PURPOSE. See the GNU Lesser General Public License or the
22 // GNU General Public License for more details.
23 //
24 // You should have received a copy of the GNU Lesser General Public
25 // License and a copy of the GNU General Public License along with
26 // Eigen. If not, see <http://www.gnu.org/licenses/>.
27 
28 #ifndef EIGEN_NUMERICAL_DIFF_H
29 #define EIGEN_NUMERICAL_DIFF_H
30 
31 namespace Eigen {
32 
33 enum NumericalDiffMode {
34  Forward,
35  Central
36 };
37 
38 
50 template<typename _Functor, NumericalDiffMode mode=Forward>
51 class NumericalDiff : public _Functor
52 {
53 public:
54  typedef _Functor Functor;
55  typedef typename Functor::Scalar Scalar;
56  typedef typename Functor::InputType InputType;
57  typedef typename Functor::ValueType ValueType;
58  typedef typename Functor::JacobianType JacobianType;
59 
60  NumericalDiff(Scalar _epsfcn=0.) : Functor(), epsfcn(_epsfcn) {}
61  NumericalDiff(const Functor& f, Scalar _epsfcn=0.) : Functor(f), epsfcn(_epsfcn) {}
62 
63  // forward constructors
64  template<typename T0>
65  NumericalDiff(const T0& a0) : Functor(a0), epsfcn(0) {}
66  template<typename T0, typename T1>
67  NumericalDiff(const T0& a0, const T1& a1) : Functor(a0, a1), epsfcn(0) {}
68  template<typename T0, typename T1, typename T2>
69  NumericalDiff(const T0& a0, const T1& a1, const T1& a2) : Functor(a0, a1, a2), epsfcn(0) {}
70 
71  enum {
72  InputsAtCompileTime = Functor::InputsAtCompileTime,
73  ValuesAtCompileTime = Functor::ValuesAtCompileTime
74  };
75 
79  int df(const InputType& _x, JacobianType &jac) const
80  {
81  /* Local variables */
82  Scalar h;
83  int nfev=0;
84  const typename InputType::Index n = _x.size();
85  const Scalar eps = internal::sqrt(((std::max)(epsfcn,NumTraits<Scalar>::epsilon() )));
86  ValueType val1, val2;
87  InputType x = _x;
88  // TODO : we should do this only if the size is not already known
89  val1.resize(Functor::values());
90  val2.resize(Functor::values());
91 
92  // initialization
93  switch(mode) {
94  case Forward:
95  // compute f(x)
96  Functor::operator()(x, val1); nfev++;
97  break;
98  case Central:
99  // do nothing
100  break;
101  default:
102  assert(false);
103  };
104 
105  // Function Body
106  for (int j = 0; j < n; ++j) {
107  h = eps * internal::abs(x[j]);
108  if (h == 0.) {
109  h = eps;
110  }
111  switch(mode) {
112  case Forward:
113  x[j] += h;
114  Functor::operator()(x, val2);
115  nfev++;
116  x[j] = _x[j];
117  jac.col(j) = (val2-val1)/h;
118  break;
119  case Central:
120  x[j] += h;
121  Functor::operator()(x, val2); nfev++;
122  x[j] -= 2*h;
123  Functor::operator()(x, val1); nfev++;
124  x[j] = _x[j];
125  jac.col(j) = (val2-val1)/(2*h);
126  break;
127  default:
128  assert(false);
129  };
130  }
131  return nfev;
132  }
133 private:
134  Scalar epsfcn;
135 
136  NumericalDiff& operator=(const NumericalDiff&);
137 };
138 
139 } // end namespace Eigen
140 
141 //vim: ai ts=4 sts=4 et sw=4
142 #endif // EIGEN_NUMERICAL_DIFF_H
143