NumericalDiff.h
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 namespace Eigen { 
00032 
00033 enum NumericalDiffMode {
00034     Forward,
00035     Central
00036 };
00037 
00038 
00050 template<typename _Functor, NumericalDiffMode mode=Forward>
00051 class NumericalDiff : public _Functor
00052 {
00053 public:
00054     typedef _Functor Functor;
00055     typedef typename Functor::Scalar Scalar;
00056     typedef typename Functor::InputType InputType;
00057     typedef typename Functor::ValueType ValueType;
00058     typedef typename Functor::JacobianType JacobianType;
00059 
00060     NumericalDiff(Scalar _epsfcn=0.) : Functor(), epsfcn(_epsfcn) {}
00061     NumericalDiff(const Functor& f, Scalar _epsfcn=0.) : Functor(f), epsfcn(_epsfcn) {}
00062 
00063     // forward constructors
00064     template<typename T0>
00065         NumericalDiff(const T0& a0) : Functor(a0), epsfcn(0) {}
00066     template<typename T0, typename T1>
00067         NumericalDiff(const T0& a0, const T1& a1) : Functor(a0, a1), epsfcn(0) {}
00068     template<typename T0, typename T1, typename T2>
00069         NumericalDiff(const T0& a0, const T1& a1, const T1& a2) : Functor(a0, a1, a2), epsfcn(0) {}
00070 
00071     enum {
00072         InputsAtCompileTime = Functor::InputsAtCompileTime,
00073         ValuesAtCompileTime = Functor::ValuesAtCompileTime
00074     };
00075 
00079     int df(const InputType& _x, JacobianType &jac) const
00080     {
00081         /* Local variables */
00082         Scalar h;
00083         int nfev=0;
00084         const typename InputType::Index n = _x.size();
00085         const Scalar eps = internal::sqrt(((std::max)(epsfcn,NumTraits<Scalar>::epsilon() )));
00086         ValueType val1, val2;
00087         InputType x = _x;
00088         // TODO : we should do this only if the size is not already known
00089         val1.resize(Functor::values());
00090         val2.resize(Functor::values());
00091 
00092         // initialization
00093         switch(mode) {
00094             case Forward:
00095                 // compute f(x)
00096                 Functor::operator()(x, val1); nfev++;
00097                 break;
00098             case Central:
00099                 // do nothing
00100                 break;
00101             default:
00102                 assert(false);
00103         };
00104 
00105         // Function Body
00106         for (int j = 0; j < n; ++j) {
00107             h = eps * internal::abs(x[j]);
00108             if (h == 0.) {
00109                 h = eps;
00110             }
00111             switch(mode) {
00112                 case Forward:
00113                     x[j] += h;
00114                     Functor::operator()(x, val2);
00115                     nfev++;
00116                     x[j] = _x[j];
00117                     jac.col(j) = (val2-val1)/h;
00118                     break;
00119                 case Central:
00120                     x[j] += h;
00121                     Functor::operator()(x, val2); nfev++;
00122                     x[j] -= 2*h;
00123                     Functor::operator()(x, val1); nfev++;
00124                     x[j] = _x[j];
00125                     jac.col(j) = (val2-val1)/(2*h);
00126                     break;
00127                 default:
00128                     assert(false);
00129             };
00130         }
00131         return nfev;
00132     }
00133 private:
00134     Scalar epsfcn;
00135 
00136     NumericalDiff& operator=(const NumericalDiff&);
00137 };
00138 
00139 } // end namespace Eigen
00140 
00141 //vim: ai ts=4 sts=4 et sw=4
00142 #endif // EIGEN_NUMERICAL_DIFF_H
00143