autodiff.cpp
Go to the documentation of this file.
00001 // This file is part of Eigen, a lightweight C++ template library
00002 // for linear algebra.
00003 //
00004 // Copyright (C) 2009 Gael Guennebaud <g.gael@free.fr>
00005 //
00006 // Eigen is free software; you can redistribute it and/or
00007 // modify it under the terms of the GNU Lesser General Public
00008 // License as published by the Free Software Foundation; either
00009 // version 3 of the License, or (at your option) any later version.
00010 //
00011 // Alternatively, you can redistribute it and/or
00012 // modify it under the terms of the GNU General Public License as
00013 // published by the Free Software Foundation; either version 2 of
00014 // the License, or (at your option) any later version.
00015 //
00016 // Eigen is distributed in the hope that it will be useful, but WITHOUT ANY
00017 // WARRANTY; without even the implied warranty of MERCHANTABILITY or FITNESS
00018 // FOR A PARTICULAR PURPOSE. See the GNU Lesser General Public License or the
00019 // GNU General Public License for more details.
00020 //
00021 // You should have received a copy of the GNU Lesser General Public
00022 // License and a copy of the GNU General Public License along with
00023 // Eigen. If not, see <http://www.gnu.org/licenses/>.
00024 
00025 #include "main.h"
00026 #include <unsupported/Eigen/AutoDiff>
00027 
00028 template<typename Scalar>
00029 EIGEN_DONT_INLINE Scalar foo(const Scalar& x, const Scalar& y)
00030 {
00031 //   return x+std::sin(y);
00032   EIGEN_ASM_COMMENT("mybegin");
00033   return static_cast<Scalar>(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));
00034   //return x+2*y*x;//x*2 -std::pow(x,2);//(2*y/x);// - y*2;
00035   EIGEN_ASM_COMMENT("myend");
00036 }
00037 
00038 template<typename _Scalar, int NX=Dynamic, int NY=Dynamic>
00039 struct TestFunc1
00040 {
00041   typedef _Scalar Scalar;
00042   enum {
00043     InputsAtCompileTime = NX,
00044     ValuesAtCompileTime = NY
00045   };
00046   typedef Matrix<Scalar,InputsAtCompileTime,1> InputType;
00047   typedef Matrix<Scalar,ValuesAtCompileTime,1> ValueType;
00048   typedef Matrix<Scalar,ValuesAtCompileTime,InputsAtCompileTime> JacobianType;
00049 
00050   int m_inputs, m_values;
00051 
00052   TestFunc1() : m_inputs(InputsAtCompileTime), m_values(ValuesAtCompileTime) {}
00053   TestFunc1(int inputs, int values) : m_inputs(inputs), m_values(values) {}
00054 
00055   int inputs() const { return m_inputs; }
00056   int values() const { return m_values; }
00057 
00058   template<typename T>
00059   void operator() (const Matrix<T,InputsAtCompileTime,1>& x, Matrix<T,ValuesAtCompileTime,1>* _v) const
00060   {
00061     Matrix<T,ValuesAtCompileTime,1>& v = *_v;
00062 
00063     v[0] = 2 * x[0] * x[0] + x[0] * x[1];
00064     v[1] = 3 * x[1] * x[0] + 0.5 * x[1] * x[1];
00065     if(inputs()>2)
00066     {
00067       v[0] += 0.5 * x[2];
00068       v[1] += x[2];
00069     }
00070     if(values()>2)
00071     {
00072       v[2] = 3 * x[1] * x[0] * x[0];
00073     }
00074     if (inputs()>2 && values()>2)
00075       v[2] *= x[2];
00076   }
00077 
00078   void operator() (const InputType& x, ValueType* v, JacobianType* _j) const
00079   {
00080     (*this)(x, v);
00081 
00082     if(_j)
00083     {
00084       JacobianType& j = *_j;
00085 
00086       j(0,0) = 4 * x[0] + x[1];
00087       j(1,0) = 3 * x[1];
00088 
00089       j(0,1) = x[0];
00090       j(1,1) = 3 * x[0] + 2 * 0.5 * x[1];
00091 
00092       if (inputs()>2)
00093       {
00094         j(0,2) = 0.5;
00095         j(1,2) = 1;
00096       }
00097       if(values()>2)
00098       {
00099         j(2,0) = 3 * x[1] * 2 * x[0];
00100         j(2,1) = 3 * x[0] * x[0];
00101       }
00102       if (inputs()>2 && values()>2)
00103       {
00104         j(2,0) *= x[2];
00105         j(2,1) *= x[2];
00106 
00107         j(2,2) = 3 * x[1] * x[0] * x[0];
00108         j(2,2) = 3 * x[1] * x[0] * x[0];
00109       }
00110     }
00111   }
00112 };
00113 
00114 template<typename Func> void forward_jacobian(const Func& f)
00115 {
00116     typename Func::InputType x = Func::InputType::Random(f.inputs());
00117     typename Func::ValueType y(f.values()), yref(f.values());
00118     typename Func::JacobianType j(f.values(),f.inputs()), jref(f.values(),f.inputs());
00119 
00120     jref.setZero();
00121     yref.setZero();
00122     f(x,&yref,&jref);
00123 //     std::cerr << y.transpose() << "\n\n";;
00124 //     std::cerr << j << "\n\n";;
00125 
00126     j.setZero();
00127     y.setZero();
00128     AutoDiffJacobian<Func> autoj(f);
00129     autoj(x, &y, &j);
00130 //     std::cerr << y.transpose() << "\n\n";;
00131 //     std::cerr << j << "\n\n";;
00132 
00133     VERIFY_IS_APPROX(y, yref);
00134     VERIFY_IS_APPROX(j, jref);
00135 }
00136 
00137 void test_autodiff_scalar()
00138 {
00139   std::cerr << foo<float>(1,2) << "\n";
00140   typedef AutoDiffScalar<Vector2f> AD;
00141   AD ax(1,Vector2f::UnitX());
00142   AD ay(2,Vector2f::UnitY());
00143   foo<AD>(ax,ay);
00144   std::cerr << foo<AD>(ax,ay).value() << " <> "
00145             << foo<AD>(ax,ay).derivatives().transpose() << "\n\n";
00146 }
00147 
00148 void test_autodiff_jacobian()
00149 {
00150   for(int i = 0; i < g_repeat; i++) {
00151     CALL_SUBTEST(( forward_jacobian(TestFunc1<double,2,2>()) ));
00152     CALL_SUBTEST(( forward_jacobian(TestFunc1<double,2,3>()) ));
00153     CALL_SUBTEST(( forward_jacobian(TestFunc1<double,3,2>()) ));
00154     CALL_SUBTEST(( forward_jacobian(TestFunc1<double,3,3>()) ));
00155     CALL_SUBTEST(( forward_jacobian(TestFunc1<double>(3,3)) ));
00156   }
00157 }
00158 
00159 void test_autodiff()
00160 {
00161     test_autodiff_scalar();
00162     test_autodiff_jacobian();
00163 }
00164 


libicr
Author(s): Robert Krug
autogenerated on Mon Jan 6 2014 11:32:29