diff options
author | Thomas Capricelli <orzel@freehackers.org> | 2009-09-28 02:43:07 +0200 |
---|---|---|
committer | Thomas Capricelli <orzel@freehackers.org> | 2009-09-28 02:43:07 +0200 |
commit | 206b5e39723c51eaa9550db04d459ebd07604415 (patch) | |
tree | 42506c765fd4021bf1fa2f9634f5d5fcf1622c6f /unsupported/Eigen | |
parent | a4532983226e04a2b5ae9f4746ca7564cbb71063 (diff) |
starting work on a Numerical differenciation module
Diffstat (limited to 'unsupported/Eigen')
-rw-r--r-- | unsupported/Eigen/NumericalDiff | 50 | ||||
-rw-r--r-- | unsupported/Eigen/src/NumericalDiff/NumericalDiff.h | 95 |
2 files changed, 145 insertions, 0 deletions
diff --git a/unsupported/Eigen/NumericalDiff b/unsupported/Eigen/NumericalDiff new file mode 100644 index 000000000..991ce7c35 --- /dev/null +++ b/unsupported/Eigen/NumericalDiff @@ -0,0 +1,50 @@ +// This file is part of Eigen, a lightweight C++ template library +// for linear algebra. +// +// Copyright (C) 2009 Thomas Capricelli <orzel@freehackers.org> +// +// Eigen is free software; you can redistribute it and/or +// modify it under the terms of the GNU Lesser General Public +// License as published by the Free Software Foundation; either +// version 3 of the License, or (at your option) any later version. +// +// Alternatively, you can redistribute it and/or +// modify it under the terms of the GNU General Public License as +// published by the Free Software Foundation; either version 2 of +// the License, or (at your option) any later version. +// +// Eigen is distributed in the hope that it will be useful, but WITHOUT ANY +// WARRANTY; without even the implied warranty of MERCHANTABILITY or FITNESS +// FOR A PARTICULAR PURPOSE. See the GNU Lesser General Public License or the +// GNU General Public License for more details. +// +// You should have received a copy of the GNU Lesser General Public +// License and a copy of the GNU General Public License along with +// Eigen. If not, see <http://www.gnu.org/licenses/>. + +#ifndef EIGEN_NUMERICALDIFF_MODULE +#define EIGEN_NUMERICALDIFF_MODULE + +#include <Eigen/Core> + +namespace Eigen { + +/** \ingroup Unsupported_modules + * \defgroup NumericalDiff_Module Support for numerical differenciation. + * See http://en.wikipedia.org/wiki/Numerical_differentiation + * + * \code + * #include <unsupported/Eigen/NumericalDiff> + * \endcode + */ +//@{ + +} + +#include "src/NumericalDiff/NumericalDiff.h" + +namespace Eigen { +//@} +} + +#endif // EIGEN_NUMERICALDIFF_MODULE diff --git a/unsupported/Eigen/src/NumericalDiff/NumericalDiff.h b/unsupported/Eigen/src/NumericalDiff/NumericalDiff.h new file mode 100644 index 000000000..223cf9e0f --- /dev/null +++ b/unsupported/Eigen/src/NumericalDiff/NumericalDiff.h @@ -0,0 +1,95 @@ +// This file is part of Eigen, a lightweight C++ template library +// for linear algebra. +// +// Copyright (C) 2009 Thomas Capricelli <orzel@freehackers.org> +// +// Eigen is free software; you can redistribute it and/or +// modify it under the terms of the GNU Lesser General Public +// License as published by the Free Software Foundation; either +// version 3 of the License, or (at your option) any later version. +// +// Alternatively, you can redistribute it and/or +// modify it under the terms of the GNU General Public License as +// published by the Free Software Foundation; either version 2 of +// the License, or (at your option) any later version. +// +// Eigen is distributed in the hope that it will be useful, but WITHOUT ANY +// WARRANTY; without even the implied warranty of MERCHANTABILITY or FITNESS +// FOR A PARTICULAR PURPOSE. See the GNU Lesser General Public License or the +// GNU General Public License for more details. +// +// You should have received a copy of the GNU Lesser General Public +// License and a copy of the GNU General Public License along with +// Eigen. If not, see <http://www.gnu.org/licenses/>. + +#ifndef EIGEN_NUMERICAL_DIFF_H +#define EIGEN_NUMERICAL_DIFF_H + +namespace Eigen +{ + +template<typename Functor> class NumericalDiff : public Functor +{ +public: + 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 T1& 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 + { + /* Local variables */ + Scalar h; + int nfev=0; + const int n = _x.size(); + const Scalar eps = ei_sqrt((std::max(epsfcn,epsilon<Scalar>() ))); + ValueType val, fx; + InputType x = _x; + // TODO : we should do this only if the size is not already known + val.resize(Functor::values()); + fx.resize(Functor::values()); + + // compute f(x) + Functor::operator()(x, fx); + + /* Function Body */ + + for (int j = 0; j < n; ++j) { + h = eps * ei_abs(x[j]); + if (h == 0.) { + h = eps; + } + x[j] += h; + Functor::operator()(x, val); + nfev++; + x[j] = _x[j]; + jac.col(j) = (val-fx)/h; + } + return nfev; + } +private: + Scalar epsfcn; +}; + +} // namespace + +#endif // EIGEN_NUMERICAL_DIFF_H |