src/NumericalDiff/NumericalDiff.h

0001 // -*- coding: utf-8
0002 // vim: set fileencoding=utf-8
0003
0004 // This file is part of Eigen, a lightweight C++ template library
0005 // for linear algebra.
0006 //
0007 // Copyright (C) 2009 Thomas Capricelli <orzel@freehackers.org>
0008 //
0009 // This Source Code Form is subject to the terms of the Mozilla
0010 // Public License v. 2.0. If a copy of the MPL was not distributed
0011 // with this file, You can obtain one at http://mozilla.org/MPL/2.0/.
0012
0013 #ifndef EIGEN_NUMERICAL_DIFF_H
0014 #define EIGEN_NUMERICAL_DIFF_H
0015
0016 namespace Eigen {
0017
0018 enum NumericalDiffMode {
0019     Forward,
0020     Central
0021 };
0022
0023
0024 /**
0025   * This class allows you to add a method df() to your functor, which will
0026   * use numerical differentiation to compute an approximate of the
0027   * derivative for the functor. Of course, if you have an analytical form
0028   * for the derivative, you should rather implement df() by yourself.
0029   *
0030   * More information on
0031   * http://en.wikipedia.org/wiki/Numerical_differentiation
0032   *
0033   * Currently only "Forward" and "Central" scheme are implemented.
0034   */
0035 template<typename _Functor, NumericalDiffMode mode=Forward>
0036 class NumericalDiff : public _Functor
0037 {
0038 public:
0039     typedef _Functor Functor;
0040     typedef typename Functor::Scalar Scalar;
0041     typedef typename Functor::InputType InputType;
0042     typedef typename Functor::ValueType ValueType;
0043     typedef typename Functor::JacobianType JacobianType;
0044
0045     NumericalDiff(Scalar _epsfcn=0.) : Functor(), epsfcn(_epsfcn) {}
0046     NumericalDiff(const Functor& f, Scalar _epsfcn=0.) : Functor(f), epsfcn(_epsfcn) {}
0047
0048     // forward constructors
0049     template<typename T0>
0050         NumericalDiff(const T0& a0) : Functor(a0), epsfcn(0) {}
0051     template<typename T0, typename T1>
0052         NumericalDiff(const T0& a0, const T1& a1) : Functor(a0, a1), epsfcn(0) {}
0053     template<typename T0, typename T1, typename T2>
0054         NumericalDiff(const T0& a0, const T1& a1, const T2& a2) : Functor(a0, a1, a2), epsfcn(0) {}
0055
0056     enum {
0057         InputsAtCompileTime = Functor::InputsAtCompileTime,
0058         ValuesAtCompileTime = Functor::ValuesAtCompileTime
0059     };
0060
0061     /**
0062       * return the number of evaluation of functor
0063      */
0064     int df(const InputType& _x, JacobianType &jac) const
0065     {
0066         using std::sqrt;
0067         using std::abs;
0068         /* Local variables */
0069         Scalar h;
0070         int nfev=0;
0071         const typename InputType::Index n = _x.size();
0072         const Scalar eps = sqrt(((std::max)(epsfcn,NumTraits<Scalar>::epsilon() )));
0073         ValueType val1, val2;
0074         InputType x = _x;
0075         // TODO : we should do this only if the size is not already known
0076         val1.resize(Functor::values());
0077         val2.resize(Functor::values());
0078
0079         // initialization
0080         switch(mode) {
0081             case Forward:
0082                 // compute f(x)
0083                 Functor::operator()(x, val1); nfev++;
0084                 break;
0085             case Central:
0086                 // do nothing
0087                 break;
0088             default:
0089                 eigen_assert(false);
0090         };
0091
0092         // Function Body
0093         for (int j = 0; j < n; ++j) {
0094             h = eps * abs(x[j]);
0095             if (h == 0.) {
0096                 h = eps;
0097             }
0098             switch(mode) {
0099                 case Forward:
0100                     x[j] += h;
0101                     Functor::operator()(x, val2);
0102                     nfev++;
0103                     x[j] = _x[j];
0104                     jac.col(j) = (val2-val1)/h;
0105                     break;
0106                 case Central:
0107                     x[j] += h;
0108                     Functor::operator()(x, val2); nfev++;
0109                     x[j] -= 2*h;
0110                     Functor::operator()(x, val1); nfev++;
0111                     x[j] = _x[j];
0112                     jac.col(j) = (val2-val1)/(2*h);
0113                     break;
0114                 default:
0115                     eigen_assert(false);
0116             };
0117         }
0118         return nfev;
0119     }
0120 private:
0121     Scalar epsfcn;
0122
0123     NumericalDiff& operator=(const NumericalDiff&);
0124 };
0125
0126 } // end namespace Eigen
0127
0128 //vim: ai ts=4 sts=4 et sw=4
0129 #endif // EIGEN_NUMERICAL_DIFF_H
0130