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 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
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
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
00089 val1.resize(Functor::values());
00090 val2.resize(Functor::values());
00091
00092
00093 switch(mode) {
00094 case Forward:
00095
00096 Functor::operator()(x, val1); nfev++;
00097 break;
00098 case Central:
00099
00100 break;
00101 default:
00102 assert(false);
00103 };
00104
00105
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 }
00140
00141
00142 #endif // EIGEN_NUMERICAL_DIFF_H
00143