| // -*- coding: utf-8 |
| // vim: set fileencoding=utf-8 |
| |
| // This file is part of Eigen, a lightweight C++ template library |
| // for linear algebra. |
| // |
| // Copyright (C) 2009 Thomas Capricelli <orzel@freehackers.org> |
| // |
| // This Source Code Form is subject to the terms of the Mozilla |
| // Public License v. 2.0. If a copy of the MPL was not distributed |
| // with this file, You can obtain one at http://mozilla.org/MPL/2.0/. |
| |
| #ifndef EIGEN_NUMERICAL_DIFF_H |
| #define EIGEN_NUMERICAL_DIFF_H |
| |
| // IWYU pragma: private |
| #include "./InternalHeaderCheck.h" |
| |
| namespace Eigen { |
| |
| enum NumericalDiffMode { Forward, Central }; |
| |
| /** |
| * This class allows you to add a method df() to your functor, which will |
| * use numerical differentiation to compute an approximate of the |
| * derivative for the functor. Of course, if you have an analytical form |
| * for the derivative, you should rather implement df() by yourself. |
| * |
| * More information on |
| * http://en.wikipedia.org/wiki/Numerical_differentiation |
| * |
| * Currently only "Forward" and "Central" scheme are implemented. |
| */ |
| template <typename Functor_, NumericalDiffMode mode = Forward> |
| class NumericalDiff : public Functor_ { |
| public: |
| typedef Functor_ Functor; |
| typedef typename Functor::Scalar Scalar; |
| typedef typename Functor::InputType InputType; |
| typedef typename Functor::ValueType ValueType; |
| typedef typename Functor::JacobianType JacobianType; |
| |
| NumericalDiff(Scalar _epsfcn = 0.) : Functor(), epsfcn(_epsfcn) {} |
| NumericalDiff(const Functor& f, Scalar _epsfcn = 0.) : Functor(f), epsfcn(_epsfcn) {} |
| |
| // forward constructors |
| template <typename T0> |
| NumericalDiff(const T0& a0) : Functor(a0), epsfcn(0) {} |
| template <typename T0, typename T1> |
| NumericalDiff(const T0& a0, const T1& a1) : Functor(a0, a1), epsfcn(0) {} |
| template <typename T0, typename T1, typename T2> |
| NumericalDiff(const T0& a0, const T1& a1, const T2& a2) : Functor(a0, a1, a2), epsfcn(0) {} |
| |
| enum { InputsAtCompileTime = Functor::InputsAtCompileTime, ValuesAtCompileTime = Functor::ValuesAtCompileTime }; |
| |
| /** |
| * return the number of evaluation of functor |
| */ |
| int df(const InputType& _x, JacobianType& jac) const { |
| using std::abs; |
| using std::sqrt; |
| /* Local variables */ |
| Scalar h; |
| int nfev = 0; |
| const typename InputType::Index n = _x.size(); |
| const Scalar eps = sqrt(((std::max)(epsfcn, NumTraits<Scalar>::epsilon()))); |
| ValueType val1, val2; |
| InputType x = _x; |
| // TODO : we should do this only if the size is not already known |
| val1.resize(Functor::values()); |
| val2.resize(Functor::values()); |
| |
| // initialization |
| switch (mode) { |
| case Forward: |
| // compute f(x) |
| Functor::operator()(x, val1); |
| nfev++; |
| break; |
| case Central: |
| // do nothing |
| break; |
| default: |
| eigen_assert(false); |
| }; |
| |
| // Function Body |
| for (int j = 0; j < n; ++j) { |
| h = eps * abs(x[j]); |
| if (h == 0.) { |
| h = eps; |
| } |
| switch (mode) { |
| case Forward: |
| x[j] += h; |
| Functor::operator()(x, val2); |
| nfev++; |
| x[j] = _x[j]; |
| jac.col(j) = (val2 - val1) / h; |
| break; |
| case Central: |
| x[j] += h; |
| Functor::operator()(x, val2); |
| nfev++; |
| x[j] -= 2 * h; |
| Functor::operator()(x, val1); |
| nfev++; |
| x[j] = _x[j]; |
| jac.col(j) = (val2 - val1) / (2 * h); |
| break; |
| default: |
| eigen_assert(false); |
| }; |
| } |
| return nfev; |
| } |
| |
| private: |
| Scalar epsfcn; |
| |
| NumericalDiff& operator=(const NumericalDiff&); |
| }; |
| |
| } // end namespace Eigen |
| |
| // vim: ai ts=4 sts=4 et sw=4 |
| #endif // EIGEN_NUMERICAL_DIFF_H |