diff options
| -rw-r--r-- | unsupported/Eigen/AutoDiff | 54 | ||||
| -rw-r--r-- | unsupported/Eigen/src/AutoDiff/AutoDiffJacobian.h | 100 | ||||
| -rw-r--r-- | unsupported/Eigen/src/AutoDiff/AutoDiffScalar.h | 321 | ||||
| -rw-r--r-- | unsupported/Eigen/src/AutoDiff/AutoDiffVector.h | 214 | ||||
| -rw-r--r-- | unsupported/Eigen/src/CMakeLists.txt | 1 | ||||
| -rw-r--r-- | unsupported/test/CMakeLists.txt | 1 | ||||
| -rw-r--r-- | unsupported/test/autodiff.cpp | 156 | ||||
| -rw-r--r-- | unsupported/test/forward_adolc.cpp | 23 |
8 files changed, 863 insertions, 7 deletions
diff --git a/unsupported/Eigen/AutoDiff b/unsupported/Eigen/AutoDiff new file mode 100644 index 000000000..f47bb0228 --- /dev/null +++ b/unsupported/Eigen/AutoDiff @@ -0,0 +1,54 @@ +// This file is part of Eigen, a lightweight C++ template library +// for linear algebra. Eigen itself is part of the KDE project. +// +// Copyright (C) 2008-2009 Gael Guennebaud <g.gael@free.fr> +// +// 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_AUTODIFF_MODULE +#define EIGEN_AUTODIFF_MODULE + +#include <Eigen/Core> + +namespace Eigen { + +/** \ingroup Unsupported_modules + * \defgroup AutoDiff_Module Auto Diff module + * + * This module features forward automatic differentation via a simple + * templated scalar type wrapper AutoDiffScalar. + * + * \code + * #include <unsupported/Eigen/AutoDiff> + * \endcode + */ +//@{ + +} + +#include "src/AutoDiff/AutoDiffScalar.h" +// #include "src/AutoDiff/AutoDiffVector.h" +#include "src/AutoDiff/AutoDiffJacobian.h" + +namespace Eigen { +//@} +} + +#endif // EIGEN_AUTODIFF_MODULE diff --git a/unsupported/Eigen/src/AutoDiff/AutoDiffJacobian.h b/unsupported/Eigen/src/AutoDiff/AutoDiffJacobian.h new file mode 100644 index 000000000..2b7987975 --- /dev/null +++ b/unsupported/Eigen/src/AutoDiff/AutoDiffJacobian.h @@ -0,0 +1,100 @@ +// This file is part of Eigen, a lightweight C++ template library +// for linear algebra. Eigen itself is part of the KDE project. +// +// Copyright (C) 2009 Gael Guennebaud <g.gael@free.fr> +// +// 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_AUTODIFF_JACOBIAN_H +#define EIGEN_AUTODIFF_JACOBIAN_H + +namespace Eigen +{ + +template<typename Functor> class AutoDiffJacobian : public Functor +{ +public: + AutoDiffJacobian() : Functor() {} + AutoDiffJacobian(const Functor& f) : Functor(f) {} + + // forward constructors + template<typename T0> + AutoDiffJacobian(const T0& a0) : Functor(a0) {} + template<typename T0, typename T1> + AutoDiffJacobian(const T0& a0, const T1& a1) : Functor(a0, a1) {} + template<typename T0, typename T1, typename T2> + AutoDiffJacobian(const T0& a0, const T1& a1, const T1& a2) : Functor(a0, a1, a2) {} + + enum { + InputsAtCompileTime = Functor::InputsAtCompileTime, + ValuesAtCompileTime = Functor::ValuesAtCompileTime + }; + + typedef typename Functor::InputType InputType; + typedef typename Functor::ValueType ValueType; + typedef typename Functor::JacobianType JacobianType; + + typedef AutoDiffScalar<Matrix<double,InputsAtCompileTime,1> > ActiveScalar; + + typedef Matrix<ActiveScalar, InputsAtCompileTime, 1> ActiveInput; + typedef Matrix<ActiveScalar, ValuesAtCompileTime, 1> ActiveValue; + + void operator() (const InputType& x, ValueType* v, JacobianType* _jac) const + { + ei_assert(v!=0); + if (!_jac) + { + Functor::operator()(x, v); + return; + } + + JacobianType& jac = *_jac; + + ActiveInput ax = x.template cast<ActiveScalar>(); + ActiveValue av(jac.rows()); + + if(InputsAtCompileTime==Dynamic) + { + for (int j=0; j<jac.cols(); j++) + ax[j].derivatives().resize(this->inputs()); + for (int j=0; j<jac.rows(); j++) + av[j].derivatives().resize(this->inputs()); + } + + for (int j=0; j<jac.cols(); j++) + for (int i=0; i<jac.cols(); i++) + ax[i].derivatives().coeffRef(j) = i==j ? 1 : 0; + + Functor::operator()(ax, &av); + + for (int i=0; i<jac.rows(); i++) + { + (*v)[i] = av[i].value(); + for (int j=0; j<jac.cols(); j++) + jac.coeffRef(i,j) = av[i].derivatives().coeff(j); + } + } +protected: + +}; + +} + +#endif // EIGEN_AUTODIFF_JACOBIAN_H diff --git a/unsupported/Eigen/src/AutoDiff/AutoDiffScalar.h b/unsupported/Eigen/src/AutoDiff/AutoDiffScalar.h new file mode 100644 index 000000000..852b43bf5 --- /dev/null +++ b/unsupported/Eigen/src/AutoDiff/AutoDiffScalar.h @@ -0,0 +1,321 @@ +// This file is part of Eigen, a lightweight C++ template library +// for linear algebra. Eigen itself is part of the KDE project. +// +// Copyright (C) 2009 Gael Guennebaud <g.gael@free.fr> +// +// 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_AUTODIFF_SCALAR_H +#define EIGEN_AUTODIFF_SCALAR_H + +namespace Eigen { + +/** \class AutoDiffScalar + * \brief A scalar type replacement with automatic differentation capability + * + * \param DerType the vector type used to store/represent the derivatives (e.g. Vector3f) + * + * This class represents a scalar value while tracking its respective derivatives. + * + * It supports the following list of global math function: + * - std::abs, std::sqrt, std::pow, std::exp, std::log, std::sin, std::cos, + * - ei_abs, ei_sqrt, ei_pow, ei_exp, ei_log, ei_sin, ei_cos, + * - ei_conj, ei_real, ei_imag, ei_abs2. + * + * AutoDiffScalar can be used as the scalar type of an Eigen::Matrix object. However, + * in that case, the expression template mechanism only occurs at the top Matrix level, + * while derivatives are computed right away. + * + */ +template<typename DerType> +class AutoDiffScalar +{ + public: + typedef typename ei_traits<DerType>::Scalar Scalar; + + inline AutoDiffScalar() {} + + inline AutoDiffScalar(const Scalar& value) + : m_value(value) + { + if(m_derivatives.size()>0) + m_derivatives.setZero(); + } + + inline AutoDiffScalar(const Scalar& value, const DerType& der) + : m_value(value), m_derivatives(der) + {} + + template<typename OtherDerType> + inline AutoDiffScalar(const AutoDiffScalar<OtherDerType>& other) + : m_value(other.value()), m_derivatives(other.derivatives()) + {} + + inline AutoDiffScalar(const AutoDiffScalar& other) + : m_value(other.value()), m_derivatives(other.derivatives()) + {} + + template<typename OtherDerType> + inline AutoDiffScalar& operator=(const AutoDiffScalar<OtherDerType>& other) + { + m_value = other.value(); + m_derivatives = other.derivatives(); + return *this; + } + + inline AutoDiffScalar& operator=(const AutoDiffScalar& other) + { + m_value = other.value(); + m_derivatives = other.derivatives(); + return *this; + } + +// inline operator const Scalar& () const { return m_value; } +// inline operator Scalar& () { return m_value; } + + inline const Scalar& value() const { return m_value; } + inline Scalar& value() { return m_value; } + + inline const DerType& derivatives() const { return m_derivatives; } + inline DerType& derivatives() { return m_derivatives; } + + template<typename OtherDerType> + inline const AutoDiffScalar<CwiseBinaryOp<ei_scalar_sum_op<Scalar>,DerType,OtherDerType> > + operator+(const AutoDiffScalar<OtherDerType>& other) const + { + return AutoDiffScalar<CwiseBinaryOp<ei_scalar_sum_op<Scalar>,DerType,OtherDerType> >( + m_value + other.value(), + m_derivatives + other.derivatives()); + } + + template<typename OtherDerType> + inline AutoDiffScalar& + operator+=(const AutoDiffScalar<OtherDerType>& other) + { + (*this) = (*this) + other; + return *this; + } + + template<typename OtherDerType> + inline const AutoDiffScalar<CwiseBinaryOp<ei_scalar_difference_op<Scalar>, DerType,OtherDerType> > + operator-(const AutoDiffScalar<OtherDerType>& other) const + { + return AutoDiffScalar<CwiseBinaryOp<ei_scalar_difference_op<Scalar>, DerType,OtherDerType> >( + m_value - other.value(), + m_derivatives - other.derivatives()); + } + + template<typename OtherDerType> + inline AutoDiffScalar& + operator-=(const AutoDiffScalar<OtherDerType>& other) + { + *this = *this - other; + return *this; + } + + template<typename OtherDerType> + inline const AutoDiffScalar<CwiseUnaryOp<ei_scalar_opposite_op<Scalar>, DerType> > + operator-() const + { + return AutoDiffScalar<CwiseUnaryOp<ei_scalar_opposite_op<Scalar>, DerType> >( + -m_value, + -m_derivatives); + } + + inline const AutoDiffScalar<CwiseUnaryOp<ei_scalar_multiple_op<Scalar>, DerType> > + operator*(const Scalar& other) const + { + return AutoDiffScalar<CwiseUnaryOp<ei_scalar_multiple_op<Scalar>, DerType> >( + m_value * other, + (m_derivatives * other)); + } + + friend inline const AutoDiffScalar<CwiseUnaryOp<ei_scalar_multiple_op<Scalar>, DerType> > + operator*(const Scalar& other, const AutoDiffScalar& a) + { + return AutoDiffScalar<CwiseUnaryOp<ei_scalar_multiple_op<Scalar>, DerType> >( + a.value() * other, + a.derivatives() * other); + } + + inline const AutoDiffScalar<CwiseUnaryOp<ei_scalar_multiple_op<Scalar>, DerType> > + operator/(const Scalar& other) const + { + return AutoDiffScalar<CwiseUnaryOp<ei_scalar_multiple_op<Scalar>, DerType> >( + m_value / other, + (m_derivatives * (Scalar(1)/other))); + } + + friend inline const AutoDiffScalar<CwiseUnaryOp<ei_scalar_multiple_op<Scalar>, DerType> > + operator/(const Scalar& other, const AutoDiffScalar& a) + { + return AutoDiffScalar<CwiseUnaryOp<ei_scalar_multiple_op<Scalar>, DerType> >( + other / a.value(), + a.derivatives() * (-Scalar(1)/other)); + } + + template<typename OtherDerType> + inline const AutoDiffScalar<CwiseUnaryOp<ei_scalar_multiple_op<Scalar>, + NestByValue<CwiseBinaryOp<ei_scalar_difference_op<Scalar>, + NestByValue<CwiseUnaryOp<ei_scalar_multiple_op<Scalar>, DerType> >, + NestByValue<CwiseUnaryOp<ei_scalar_multiple_op<Scalar>, OtherDerType> > > > > > + operator/(const AutoDiffScalar<OtherDerType>& other) const + { + return AutoDiffScalar<CwiseUnaryOp<ei_scalar_multiple_op<Scalar>, + NestByValue<CwiseBinaryOp<ei_scalar_difference_op<Scalar>, + NestByValue<CwiseUnaryOp<ei_scalar_multiple_op<Scalar>, DerType> >, + NestByValue<CwiseUnaryOp<ei_scalar_multiple_op<Scalar>, OtherDerType> > > > > >( + m_value / other.value(), + ((m_derivatives * other.value()).nestByValue() - (m_value * other.derivatives()).nestByValue()).nestByValue() + * (Scalar(1)/(other.value()*other.value()))); + } + + template<typename OtherDerType> + inline const AutoDiffScalar<CwiseBinaryOp<ei_scalar_sum_op<Scalar>, + NestByValue<CwiseUnaryOp<ei_scalar_multiple_op<Scalar>, DerType> >, + NestByValue<CwiseUnaryOp<ei_scalar_multiple_op<Scalar>, OtherDerType> > > > + operator*(const AutoDiffScalar<OtherDerType>& other) const + { + return AutoDiffScalar<CwiseBinaryOp<ei_scalar_sum_op<Scalar>, + NestByValue<CwiseUnaryOp<ei_scalar_multiple_op<Scalar>, DerType> >, + NestByValue<CwiseUnaryOp<ei_scalar_multiple_op<Scalar>, OtherDerType> > > >( + m_value * other.value(), + (m_derivatives * other.value()).nestByValue() + (m_value * other.derivatives()).nestByValue()); + } + + inline AutoDiffScalar& operator*=(const Scalar& other) + { + *this = *this * other; + return *this; + } + + template<typename OtherDerType> + inline AutoDiffScalar& operator*=(const AutoDiffScalar<OtherDerType>& other) + { + *this = *this * other; + return *this; + } + + protected: + Scalar m_value; + DerType m_derivatives; + +}; + +} + +#define EIGEN_AUTODIFF_DECLARE_GLOBAL_UNARY(FUNC,CODE) \ + template<typename DerType> \ + inline const AutoDiffScalar<CwiseUnaryOp<ei_scalar_multiple_op<typename ei_traits<DerType>::Scalar>, DerType> > \ + FUNC(const AutoDiffScalar<DerType>& x) { \ + typedef typename ei_traits<DerType>::Scalar Scalar; \ + typedef AutoDiffScalar<CwiseUnaryOp<ei_scalar_multiple_op<Scalar>, DerType> > ReturnType; \ + CODE; \ + } + +namespace std +{ + EIGEN_AUTODIFF_DECLARE_GLOBAL_UNARY(abs, + return ReturnType(std::abs(x.value()), x.derivatives() * (sign(x.value())));) + + EIGEN_AUTODIFF_DECLARE_GLOBAL_UNARY(sqrt, + Scalar sqrtx = std::sqrt(x.value()); + return ReturnType(sqrtx,x.derivatives() * (Scalar(0.5) / sqrtx));) + + EIGEN_AUTODIFF_DECLARE_GLOBAL_UNARY(cos, + return ReturnType(std::cos(x.value()), x.derivatives() * (-std::sin(x.value())));) + + EIGEN_AUTODIFF_DECLARE_GLOBAL_UNARY(sin, + return ReturnType(std::sin(x.value()),x.derivatives() * std::cos(x.value()));) + + EIGEN_AUTODIFF_DECLARE_GLOBAL_UNARY(exp, + Scalar expx = std::exp(x.value()); + return ReturnType(expx,x.derivatives() * expx);) + + EIGEN_AUTODIFF_DECLARE_GLOBAL_UNARY(ei_log, + return ReturnType(std::log(x.value),x.derivatives() * (Scalar(1).x.value()));) + + template<typename DerType> + inline const AutoDiffScalar<CwiseUnaryOp<ei_scalar_multiple_op<typename ei_traits<DerType>::Scalar>, DerType> > + pow(const AutoDiffScalar<DerType>& x, typename ei_traits<DerType>::Scalar y) + { + typedef typename ei_traits<DerType>::Scalar Scalar; + return AutoDiffScalar<CwiseUnaryOp<ei_scalar_multiple_op<Scalar>, DerType> >( + std::pow(x.value(),y), + x.derivatives() * (y * std::pow(x.value(),y-1))); + } + +} + +namespace Eigen { + +template<typename DerType> +inline const AutoDiffScalar<DerType>& ei_conj(const AutoDiffScalar<DerType>& x) { return x; } +template<typename DerType> +inline const AutoDiffScalar<DerType>& ei_real(const AutoDiffScalar<DerType>& x) { return x; } +template<typename DerType> +inline typename DerType::Scalar ei_imag(const AutoDiffScalar<DerType>&) { return 0.; } + +EIGEN_AUTODIFF_DECLARE_GLOBAL_UNARY(ei_abs, + return ReturnType(ei_abs(x.value()), x.derivatives() * (sign(x.value())));) + +EIGEN_AUTODIFF_DECLARE_GLOBAL_UNARY(ei_abs2, + return ReturnType(ei_abs2(x.value()), x.derivatives() * (Scalar(2)*x.value()));) + +EIGEN_AUTODIFF_DECLARE_GLOBAL_UNARY(ei_sqrt, + Scalar sqrtx = ei_sqrt(x.value()); + return ReturnType(sqrtx,x.derivatives() * (Scalar(0.5) / sqrtx));) + +EIGEN_AUTODIFF_DECLARE_GLOBAL_UNARY(ei_cos, + return ReturnType(ei_cos(x.value()), x.derivatives() * (-ei_sin(x.value())));) + +EIGEN_AUTODIFF_DECLARE_GLOBAL_UNARY(ei_sin, + return ReturnType(ei_sin(x.value()),x.derivatives() * ei_cos(x.value()));) + +EIGEN_AUTODIFF_DECLARE_GLOBAL_UNARY(ei_exp, + Scalar expx = ei_exp(x.value()); + return ReturnType(expx,x.derivatives() * expx);) + +EIGEN_AUTODIFF_DECLARE_GLOBAL_UNARY(ei_log, + return ReturnType(ei_log(x.value),x.derivatives() * (Scalar(1).x.value()));) + +template<typename DerType> +inline const AutoDiffScalar<CwiseUnaryOp<ei_scalar_multiple_op<typename ei_traits<DerType>::Scalar>, DerType> > +ei_pow(const AutoDiffScalar<DerType>& x, typename ei_traits<DerType>::Scalar y) +{ return std::pow(x,y);} + +#undef EIGEN_AUTODIFF_DECLARE_GLOBAL_UNARY + +template<typename DerType> struct NumTraits<AutoDiffScalar<DerType> > +{ + typedef typename DerType::Scalar Real; + typedef AutoDiffScalar<DerType> FloatingPoint; + enum { + IsComplex = 0, + HasFloatingPoint = 1, + ReadCost = 1, + AddCost = 1, + MulCost = 1 + }; +}; + +} + +#endif // EIGEN_AUTODIFF_SCALAR_H diff --git a/unsupported/Eigen/src/AutoDiff/AutoDiffVector.h b/unsupported/Eigen/src/AutoDiff/AutoDiffVector.h new file mode 100644 index 000000000..a58bb6b2e --- /dev/null +++ b/unsupported/Eigen/src/AutoDiff/AutoDiffVector.h @@ -0,0 +1,214 @@ +// This file is part of Eigen, a lightweight C++ template library +// for linear algebra. Eigen itself is part of the KDE project. +// +// Copyright (C) 2009 Gael Guennebaud <g.gael@free.fr> +// +// 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_AUTODIFF_VECTOR_H +#define EIGEN_AUTODIFF_VECTOR_H + +namespace Eigen { + +/* \class AutoDiffScalar + * \brief A scalar type replacement with automatic differentation capability + * + * \param DerType the vector type used to store/represent the derivatives (e.g. Vector3f) + * + * This class represents a scalar value while tracking its respective derivatives. + * + * It supports the following list of global math function: + * - std::abs, std::sqrt, std::pow, std::exp, std::log, std::sin, std::cos, + * - ei_abs, ei_sqrt, ei_pow, ei_exp, ei_log, ei_sin, ei_cos, + * - ei_conj, ei_real, ei_imag, ei_abs2. + * + * AutoDiffScalar can be used as the scalar type of an Eigen::Matrix object. However, + * in that case, the expression template mechanism only occurs at the top Matrix level, + * while derivatives are computed right away. + * + */ +template<typename ValueType, typename JacobianType> +class AutoDiffVector +{ + public: + typedef typename ei_traits<ValueType>::Scalar Scalar; + + inline AutoDiffVector() {} + + inline AutoDiffVector(const ValueType& values) + : m_values(values) + { + m_jacobian.setZero(); + } + + inline AutoDiffVector(const ValueType& values, const JacobianType& jac) + : m_values(values), m_jacobian(jac) + {} + + template<typename OtherValueType, typename OtherJacobianType> + inline AutoDiffVector(const AutoDiffVector<OtherValueType, OtherJacobianType>& other) + : m_values(other.values()), m_jacobian(other.jacobian()) + {} + + inline AutoDiffVector(const AutoDiffVector& other) + : m_values(other.values()), m_jacobian(other.jacobian()) + {} + + template<typename OtherValueType, typename OtherJacobianType> + inline AutoDiffScalar& operator=(const AutoDiffVector<OtherValueType, OtherJacobianType>& other) + { + m_values = other.values(); + m_jacobian = other.jacobian(); + return *this; + } + + inline AutoDiffVector& operator=(const AutoDiffVector& other) + { + m_values = other.values(); + m_jacobian = other.jacobian(); + return *this; + } + + inline const ValueType& values() const { return m_values; } + inline ValueType& values() { return m_values; } + + inline const JacobianType& jacobian() const { return m_jacobian; } + inline JacobianType& jacobian() { return m_jacobian; } + + template<typename OtherValueType,typename OtherJacobianType> + inline const AutoDiffVector< + CwiseBinaryOp<ei_scalar_sum_op<Scalar>,ValueType,OtherValueType> > + CwiseBinaryOp<ei_scalar_sum_op<Scalar>,JacobianType,OtherJacobianType> > + operator+(const AutoDiffScalar<OtherDerType>& other) const + { + return AutoDiffVector< + CwiseBinaryOp<ei_scalar_sum_op<Scalar>,ValueType,OtherValueType> > + CwiseBinaryOp<ei_scalar_sum_op<Scalar>,JacobianType,OtherJacobianType> >( + m_values + other.values(), + m_jacobian + other.jacobian()); + } + + template<typename OtherValueType, typename OtherJacobianType> + inline AutoDiffVector& + operator+=(const AutoDiffVector<OtherValueType,OtherDerType>& other) + { + m_values += other.values(); + m_jacobian += other.jacobian(); + return *this; + } + + template<typename OtherValueType,typename OtherJacobianType> + inline const AutoDiffVector< + CwiseBinaryOp<ei_scalar_difference_op<Scalar>,ValueType,OtherValueType> > + CwiseBinaryOp<ei_scalar_difference_op<Scalar>,JacobianType,OtherJacobianType> > + operator-(const AutoDiffScalar<OtherDerType>& other) const + { + return AutoDiffVector< + CwiseBinaryOp<ei_scalar_difference_op<Scalar>,ValueType,OtherValueType> > + CwiseBinaryOp<ei_scalar_difference_op<Scalar>,JacobianType,OtherJacobianType> >( + m_values - other.values(), + m_jacobian - other.jacobian()); + } + + template<typename OtherValueType, typename OtherJacobianType> + inline AutoDiffVector& + operator-=(const AutoDiffVector<OtherValueType,OtherDerType>& other) + { + m_values -= other.values(); + m_jacobian -= other.jacobian(); + return *this; + } + + inline const AutoDiffVector< + CwiseUnaryOp<ei_scalar_opposite_op<Scalar>, ValueType> + CwiseUnaryOp<ei_scalar_opposite_op<Scalar>, JacobianType> > + operator-() const + { + return AutoDiffVector< + CwiseUnaryOp<ei_scalar_opposite_op<Scalar>, ValueType> + CwiseUnaryOp<ei_scalar_opposite_op<Scalar>, JacobianType> >( + -m_values, + -m_jacobian); + } + + inline const AutoDiffVector< + CwiseUnaryOp<ei_scalar_multiple_op<Scalar>, ValueType> + CwiseUnaryOp<ei_scalar_multiple_op<Scalar>, JacobianType> > + operator*(const Scalar& other) const + { + return AutoDiffVector< + CwiseUnaryOp<ei_scalar_multiple_op<Scalar>, ValueType> + CwiseUnaryOp<ei_scalar_multiple_op<Scalar>, JacobianType> >( + m_values * other, + (m_jacobian * other)); + } + + friend inline const AutoDiffVector< + CwiseUnaryOp<ei_scalar_multiple_op<Scalar>, ValueType> + CwiseUnaryOp<ei_scalar_multiple_op<Scalar>, JacobianType> > + operator*(const Scalar& other, const AutoDiffVector& v) + { + return AutoDiffVector< + CwiseUnaryOp<ei_scalar_multiple_op<Scalar>, ValueType> + CwiseUnaryOp<ei_scalar_multiple_op<Scalar>, JacobianType> >( + v.values() * other, + v.jacobian() * other); + } + +// template<typename OtherValueType,typename OtherJacobianType> +// inline const AutoDiffVector< +// CwiseBinaryOp<ei_scalar_multiple_op<Scalar>, ValueType, OtherValueType> +// CwiseBinaryOp<ei_scalar_sum_op<Scalar>, +// NestByValue<CwiseUnaryOp<ei_scalar_multiple_op<Scalar>, JacobianType> >, +// NestByValue<CwiseUnaryOp<ei_scalar_multiple_op<Scalar>, OtherJacobianType> > > > +// operator*(const AutoDiffVector<OtherValueType,OtherJacobianType>& other) const +// { +// return AutoDiffVector< +// CwiseBinaryOp<ei_scalar_multiple_op<Scalar>, ValueType, OtherValueType> +// CwiseBinaryOp<ei_scalar_sum_op<Scalar>, +// NestByValue<CwiseUnaryOp<ei_scalar_multiple_op<Scalar>, JacobianType> >, +// NestByValue<CwiseUnaryOp<ei_scalar_multiple_op<Scalar>, OtherJacobianType> > > >( +// m_values.cwise() * other.values(), +// (m_jacobian * other.values()).nestByValue() + (m_values * other.jacobian()).nestByValue()); +// } + + inline AutoDiffVector& operator*=(const Scalar& other) + { + m_values *= other; + m_jacobian *= other; + return *this; + } + + template<typename OtherValueType,typename OtherJacobianType> + inline AutoDiffVector& operator*=(const AutoDiffVector<OtherValueType,OtherJacobianType>& other) + { + *this = *this * other; + return *this; + } + + protected: + ValueType m_values; + JacobianType m_jacobian; + +}; + +} + +#endif // EIGEN_AUTODIFF_VECTOR_H diff --git a/unsupported/Eigen/src/CMakeLists.txt b/unsupported/Eigen/src/CMakeLists.txt index 727b18cf5..c2f63db20 100644 --- a/unsupported/Eigen/src/CMakeLists.txt +++ b/unsupported/Eigen/src/CMakeLists.txt @@ -1,2 +1,3 @@ ADD_SUBDIRECTORY(IterativeSolvers) ADD_SUBDIRECTORY(BVH) +ADD_SUBDIRECTORY(AutoDiff) diff --git a/unsupported/test/CMakeLists.txt b/unsupported/test/CMakeLists.txt index 2e4bde3e1..16ad1e7a1 100644 --- a/unsupported/test/CMakeLists.txt +++ b/unsupported/test/CMakeLists.txt @@ -15,4 +15,5 @@ else(ADOLC_FOUND) ei_add_property(EIGEN_MISSING_BACKENDS "Adolc") endif(ADOLC_FOUND) +ei_add_test(autodiff) ei_add_test(BVH) diff --git a/unsupported/test/autodiff.cpp b/unsupported/test/autodiff.cpp new file mode 100644 index 000000000..7d619897c --- /dev/null +++ b/unsupported/test/autodiff.cpp @@ -0,0 +1,156 @@ +// This file is part of Eigen, a lightweight C++ template library +// for linear algebra. Eigen itself is part of the KDE project. +// +// Copyright (C) 2009 Gael Guennebaud <g.gael@free.fr> +// +// 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/>. + +#include "main.h" +#include <unsupported/Eigen/AutoDiff> + +template<typename Scalar> +EIGEN_DONT_INLINE Scalar foo(const Scalar& x, const Scalar& y) +{ +// return x+std::sin(y); + asm("#mybegin"); + return x*2 - std::pow(x,2) + 2*std::sqrt(y*y) - 4 * std::sin(x) + 2 * std::cos(y) - std::exp(-0.5*x*x); +// return y/x;// - y*2; + asm("#myend"); +} + +template<typename _Scalar, int NX=Dynamic, int NY=Dynamic> +struct TestFunc1 +{ + typedef _Scalar Scalar; + enum { + InputsAtCompileTime = NX, + ValuesAtCompileTime = NY + }; + typedef Matrix<Scalar,InputsAtCompileTime,1> InputType; + typedef Matrix<Scalar,ValuesAtCompileTime,1> ValueType; + typedef Matrix<Scalar,ValuesAtCompileTime,InputsAtCompileTime> JacobianType; + + int m_inputs, m_values; + + TestFunc1() : m_inputs(InputsAtCompileTime), m_values(ValuesAtCompileTime) {} + TestFunc1(int inputs, int values) : m_inputs(inputs), m_values(values) {} + + int inputs() const { return m_inputs; } + int values() const { return m_values; } + + template<typename T> + void operator() (const Matrix<T,InputsAtCompileTime,1>& x, Matrix<T,ValuesAtCompileTime,1>* _v) const + { + Matrix<T,ValuesAtCompileTime,1>& v = *_v; + + v[0] = 2 * x[0] * x[0] + x[0] * x[1]; + v[1] = 3 * x[1] * x[0] + 0.5 * x[1] * x[1]; + if(inputs()>2) + { + v[0] += 0.5 * x[2]; + v[1] += x[2]; + } + if(values()>2) + { + v[2] = 3 * x[1] * x[0] * x[0]; + } + if (inputs()>2 && values()>2) + v[2] *= x[2]; + } + + void operator() (const InputType& x, ValueType* v, JacobianType* _j) const + { + (*this)(x, v); + + if(_j) + { + JacobianType& j = *_j; + + j(0,0) = 4 * x[0] + x[1]; + j(1,0) = 3 * x[1]; + + j(0,1) = x[0]; + j(1,1) = 3 * x[0] + 2 * 0.5 * x[1]; + + if (inputs()>2) + { + j(0,2) = 0.5; + j(1,2) = 1; + } + if(values()>2) + { + j(2,0) = 3 * x[1] * 2 * x[0]; + j(2,1) = 3 * x[0] * x[0]; + } + if (inputs()>2 && values()>2) + { + j(2,0) *= x[2]; + j(2,1) *= x[2]; + + j(2,2) = 3 * x[1] * x[0] * x[0]; + j(2,2) = 3 * x[1] * x[0] * x[0]; + } + } + } +}; + +template<typename Func> void adolc_forward_jacobian(const Func& f) +{ + typename Func::InputType x = Func::InputType::Random(f.inputs()); + typename Func::ValueType y(f.values()), yref(f.values()); + typename Func::JacobianType j(f.values(),f.inputs()), jref(f.values(),f.inputs()); + + jref.setZero(); + yref.setZero(); + f(x,&yref,&jref); +// std::cerr << y.transpose() << "\n\n";; +// std::cerr << j << "\n\n";; + + j.setZero(); + y.setZero(); + AutoDiffJacobian<Func> autoj(f); + autoj(x, &y, &j); +// std::cerr << y.transpose() << "\n\n";; +// std::cerr << j << "\n\n";; + + VERIFY_IS_APPROX(y, yref); + VERIFY_IS_APPROX(j, jref); +} + +void test_autodiff() +{ + std::sqrt(3); + std::sin(3); + std::cerr << foo<float>(1,2) << "\n"; + AutoDiffScalar<Vector2f> ax(1,Vector2f::UnitX()); + AutoDiffScalar<Vector2f> ay(2,Vector2f::UnitY()); + std::cerr << foo<AutoDiffScalar<Vector2f> >(ax,ay).value() << " <> " + << foo<AutoDiffScalar<Vector2f> >(ax,ay).derivatives().transpose() << "\n\n"; + + for(int i = 0; i < g_repeat; i++) { + CALL_SUBTEST(( adolc_forward_jacobian(TestFunc1<double,2,2>()) )); + CALL_SUBTEST(( adolc_forward_jacobian(TestFunc1<double,2,3>()) )); + CALL_SUBTEST(( adolc_forward_jacobian(TestFunc1<double,3,2>()) )); + CALL_SUBTEST(( adolc_forward_jacobian(TestFunc1<double,3,3>()) )); + CALL_SUBTEST(( adolc_forward_jacobian(TestFunc1<double>(3,3)) )); + } + +// exit(1); +} diff --git a/unsupported/test/forward_adolc.cpp b/unsupported/test/forward_adolc.cpp index 016e20cdb..aa7618bff 100644 --- a/unsupported/test/forward_adolc.cpp +++ b/unsupported/test/forward_adolc.cpp @@ -28,7 +28,7 @@ int adtl::ADOLC_numDir; -template<typename _Scalar, int NX, int NY> +template<typename _Scalar, int NX=Dynamic, int NY=Dynamic> struct TestFunc1 { typedef _Scalar Scalar; @@ -39,6 +39,14 @@ struct TestFunc1 typedef Matrix<Scalar,InputsAtCompileTime,1> InputType; typedef Matrix<Scalar,ValuesAtCompileTime,1> ValueType; typedef Matrix<Scalar,ValuesAtCompileTime,InputsAtCompileTime> JacobianType; + + int m_inputs, m_values; + + TestFunc1() : m_inputs(InputsAtCompileTime), m_values(ValuesAtCompileTime) {} + TestFunc1(int inputs, int values) : m_inputs(inputs), m_values(values) {} + + int inputs() const { return m_inputs; } + int values() const { return m_values } template<typename T> void operator() (const Matrix<T,InputsAtCompileTime,1>& x, Matrix<T,ValuesAtCompileTime,1>* _v) const @@ -47,16 +55,16 @@ struct TestFunc1 v[0] = 2 * x[0] * x[0] + x[0] * x[1]; v[1] = 3 * x[1] * x[0] + 0.5 * x[1] * x[1]; - if(NX>2) + if(inputs()>2) { v[0] += 0.5 * x[2]; v[1] += x[2]; } - if(NY>2) + if(values()>2) { v[2] = 3 * x[1] * x[0] * x[0]; } - if (NX>2 && NY>2) + if (inputs()>2 && values()>2) v[2] *= x[2]; } @@ -74,17 +82,17 @@ struct TestFunc1 j(0,1) = x[0]; j(1,1) = 3 * x[0] + 2 * 0.5 * x[1]; - if (NX>2) + if (inputs()>2) { j(0,2) = 0.5; j(1,2) = 1; } - if(NY>2) + if(values()>2) { j(2,0) = 3 * x[1] * 2 * x[0]; j(2,1) = 3 * x[0] * x[0]; } - if (NX>2 && NY>2) + if (inputs()>2 && values()>2) { j(2,0) *= x[2]; j(2,1) *= x[2]; @@ -128,5 +136,6 @@ void test_forward_adolc() CALL_SUBTEST(( adolc_forward_jacobian(TestFunc1<double,2,3>()) )); CALL_SUBTEST(( adolc_forward_jacobian(TestFunc1<double,3,2>()) )); CALL_SUBTEST(( adolc_forward_jacobian(TestFunc1<double,3,3>()) )); + CALL_SUBTEST(( adolc_forward_jacobian(TestFunc1<double>(3,3)) )); } } |