NumericalDiff.h
Go to the documentation of this file.
00001 // -*- coding: utf-8
00002 // vim: set fileencoding=utf-8
00003 
00004 // This file is part of Eigen, a lightweight C++ template library
00005 // for linear algebra.
00006 //
00007 // Copyright (C) 2009 Thomas Capricelli <orzel@freehackers.org>
00008 //
00009 // Eigen is free software; you can redistribute it and/or
00010 // modify it under the terms of the GNU Lesser General Public
00011 // License as published by the Free Software Foundation; either
00012 // version 3 of the License, or (at your option) any later version.
00013 //
00014 // Alternatively, you can redistribute it and/or
00015 // modify it under the terms of the GNU General Public License as
00016 // published by the Free Software Foundation; either version 2 of
00017 // the License, or (at your option) any later version.
00018 //
00019 // Eigen is distributed in the hope that it will be useful, but WITHOUT ANY
00020 // WARRANTY; without even the implied warranty of MERCHANTABILITY or FITNESS
00021 // FOR A PARTICULAR PURPOSE. See the GNU Lesser General Public License or the
00022 // GNU General Public License for more details.
00023 //
00024 // You should have received a copy of the GNU Lesser General Public
00025 // License and a copy of the GNU General Public License along with
00026 // Eigen. If not, see <http://www.gnu.org/licenses/>.
00027 
00028 #ifndef EIGEN_NUMERICAL_DIFF_H
00029 #define EIGEN_NUMERICAL_DIFF_H
00030 
00031 enum NumericalDiffMode {
00032     Forward,
00033     Central
00034 };
00035 
00036 
00048 template<typename _Functor, NumericalDiffMode mode=Forward>
00049 class NumericalDiff : public _Functor
00050 {
00051 public:
00052     typedef _Functor Functor;
00053     typedef typename Functor::Scalar Scalar;
00054     typedef typename Functor::InputType InputType;
00055     typedef typename Functor::ValueType ValueType;
00056     typedef typename Functor::JacobianType JacobianType;
00057 
00058     NumericalDiff(Scalar _epsfcn=0.) : Functor(), epsfcn(_epsfcn) {}
00059     NumericalDiff(const Functor& f, Scalar _epsfcn=0.) : Functor(f), epsfcn(_epsfcn) {}
00060 
00061     // forward constructors
00062     template<typename T0>
00063         NumericalDiff(const T0& a0) : Functor(a0), epsfcn(0) {}
00064     template<typename T0, typename T1>
00065         NumericalDiff(const T0& a0, const T1& a1) : Functor(a0, a1), epsfcn(0) {}
00066     template<typename T0, typename T1, typename T2>
00067         NumericalDiff(const T0& a0, const T1& a1, const T1& a2) : Functor(a0, a1, a2), epsfcn(0) {}
00068 
00069     enum {
00070         InputsAtCompileTime = Functor::InputsAtCompileTime,
00071         ValuesAtCompileTime = Functor::ValuesAtCompileTime
00072     };
00073 
00077     int df(const InputType& _x, JacobianType &jac) const
00078     {
00079         /* Local variables */
00080         Scalar h;
00081         int nfev=0;
00082         const typename InputType::Index n = _x.size();
00083         const Scalar eps = internal::sqrt((std::max(epsfcn,NumTraits<Scalar>::epsilon() )));
00084         ValueType val1, val2;
00085         InputType x = _x;
00086         // TODO : we should do this only if the size is not already known
00087         val1.resize(Functor::values());
00088         val2.resize(Functor::values());
00089 
00090         // initialization
00091         switch(mode) {
00092             case Forward:
00093                 // compute f(x)
00094                 Functor::operator()(x, val1); nfev++;
00095                 break;
00096             case Central:
00097                 // do nothing
00098                 break;
00099             default:
00100                 assert(false);
00101         };
00102 
00103         // Function Body
00104         for (int j = 0; j < n; ++j) {
00105             h = eps * internal::abs(x[j]);
00106             if (h == 0.) {
00107                 h = eps;
00108             }
00109             switch(mode) {
00110                 case Forward:
00111                     x[j] += h;
00112                     Functor::operator()(x, val2);
00113                     nfev++;
00114                     x[j] = _x[j];
00115                     jac.col(j) = (val2-val1)/h;
00116                     break;
00117                 case Central:
00118                     x[j] += h;
00119                     Functor::operator()(x, val2); nfev++;
00120                     x[j] -= 2*h;
00121                     Functor::operator()(x, val1); nfev++;
00122                     x[j] = _x[j];
00123                     jac.col(j) = (val2-val1)/(2*h);
00124                     break;
00125                 default:
00126                     assert(false);
00127             };
00128         }
00129         return nfev;
00130     }
00131 private:
00132     Scalar epsfcn;
00133 
00134     NumericalDiff& operator=(const NumericalDiff&);
00135 };
00136 
00137 //vim: ai ts=4 sts=4 et sw=4
00138 #endif // EIGEN_NUMERICAL_DIFF_H
00139 


re_vision
Author(s): Dorian Galvez-Lopez
autogenerated on Sun Jan 5 2014 11:32:01