Go to the documentation of this file.00001
00002
00003
00004
00005
00006
00007
00008
00009
00010
00011
00012
00013
00014
00015
00016
00017
00018
00019
00020
00021
00022
00023
00024
00025
00026
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
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
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
00087 val1.resize(Functor::values());
00088 val2.resize(Functor::values());
00089
00090
00091 switch(mode) {
00092 case Forward:
00093
00094 Functor::operator()(x, val1); nfev++;
00095 break;
00096 case Central:
00097
00098 break;
00099 default:
00100 assert(false);
00101 };
00102
00103
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
00138 #endif // EIGEN_NUMERICAL_DIFF_H
00139