1c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath// -*- coding: utf-8 2c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath// vim: set fileencoding=utf-8 3c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath 4c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath// This file is part of Eigen, a lightweight C++ template library 5c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath// for linear algebra. 6c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath// 7c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath// Copyright (C) 2009 Thomas Capricelli <orzel@freehackers.org> 8c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath// 9c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath// This Source Code Form is subject to the terms of the Mozilla 10c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath// Public License v. 2.0. If a copy of the MPL was not distributed 11c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath// with this file, You can obtain one at http://mozilla.org/MPL/2.0/. 12c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath 13c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath#ifndef EIGEN_NUMERICAL_DIFF_H 14c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath#define EIGEN_NUMERICAL_DIFF_H 15c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath 16c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamathnamespace Eigen { 17c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath 18c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamathenum NumericalDiffMode { 19c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath Forward, 20c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath Central 21c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath}; 22c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath 23c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath 24c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath/** 25c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath * This class allows you to add a method df() to your functor, which will 26c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath * use numerical differentiation to compute an approximate of the 27c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath * derivative for the functor. Of course, if you have an analytical form 28c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath * for the derivative, you should rather implement df() by yourself. 29c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath * 30c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath * More information on 31c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath * http://en.wikipedia.org/wiki/Numerical_differentiation 32c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath * 33c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath * Currently only "Forward" and "Central" scheme are implemented. 34c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath */ 35c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamathtemplate<typename _Functor, NumericalDiffMode mode=Forward> 36c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamathclass NumericalDiff : public _Functor 37c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath{ 38c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamathpublic: 39c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath typedef _Functor Functor; 40c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath typedef typename Functor::Scalar Scalar; 41c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath typedef typename Functor::InputType InputType; 42c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath typedef typename Functor::ValueType ValueType; 43c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath typedef typename Functor::JacobianType JacobianType; 44c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath 45c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath NumericalDiff(Scalar _epsfcn=0.) : Functor(), epsfcn(_epsfcn) {} 46c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath NumericalDiff(const Functor& f, Scalar _epsfcn=0.) : Functor(f), epsfcn(_epsfcn) {} 47c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath 48c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath // forward constructors 49c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath template<typename T0> 50c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath NumericalDiff(const T0& a0) : Functor(a0), epsfcn(0) {} 51c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath template<typename T0, typename T1> 52c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath NumericalDiff(const T0& a0, const T1& a1) : Functor(a0, a1), epsfcn(0) {} 53c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath template<typename T0, typename T1, typename T2> 54c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath NumericalDiff(const T0& a0, const T1& a1, const T2& a2) : Functor(a0, a1, a2), epsfcn(0) {} 55c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath 56c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath enum { 57c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath InputsAtCompileTime = Functor::InputsAtCompileTime, 58c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath ValuesAtCompileTime = Functor::ValuesAtCompileTime 59c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath }; 60c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath 61c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath /** 62c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath * return the number of evaluation of functor 63c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath */ 64c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath int df(const InputType& _x, JacobianType &jac) const 65c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath { 667faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez using std::sqrt; 677faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez using std::abs; 68c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath /* Local variables */ 69c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath Scalar h; 70c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath int nfev=0; 71c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath const typename InputType::Index n = _x.size(); 727faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez const Scalar eps = sqrt(((std::max)(epsfcn,NumTraits<Scalar>::epsilon() ))); 73c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath ValueType val1, val2; 74c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath InputType x = _x; 75c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath // TODO : we should do this only if the size is not already known 76c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath val1.resize(Functor::values()); 77c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath val2.resize(Functor::values()); 78c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath 79c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath // initialization 80c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath switch(mode) { 81c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath case Forward: 82c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath // compute f(x) 83c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath Functor::operator()(x, val1); nfev++; 84c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath break; 85c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath case Central: 86c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath // do nothing 87c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath break; 88c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath default: 897faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez eigen_assert(false); 90c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath }; 91c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath 92c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath // Function Body 93c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath for (int j = 0; j < n; ++j) { 947faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez h = eps * abs(x[j]); 95c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath if (h == 0.) { 96c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath h = eps; 97c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath } 98c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath switch(mode) { 99c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath case Forward: 100c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath x[j] += h; 101c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath Functor::operator()(x, val2); 102c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath nfev++; 103c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath x[j] = _x[j]; 104c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath jac.col(j) = (val2-val1)/h; 105c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath break; 106c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath case Central: 107c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath x[j] += h; 108c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath Functor::operator()(x, val2); nfev++; 109c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath x[j] -= 2*h; 110c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath Functor::operator()(x, val1); nfev++; 111c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath x[j] = _x[j]; 112c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath jac.col(j) = (val2-val1)/(2*h); 113c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath break; 114c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath default: 1157faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez eigen_assert(false); 116c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath }; 117c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath } 118c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath return nfev; 119c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath } 120c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamathprivate: 121c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath Scalar epsfcn; 122c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath 123c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath NumericalDiff& operator=(const NumericalDiff&); 124c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath}; 125c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath 126c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath} // end namespace Eigen 127c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath 128c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath//vim: ai ts=4 sts=4 et sw=4 129c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath#endif // EIGEN_NUMERICAL_DIFF_H 130c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath 131