forward_adolc.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) 2008 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 #define NUMBER_DIRECTIONS 16
00027 #include <unsupported/Eigen/AdolcForward>
00028 
00029 int adtl::ADOLC_numDir;
00030 
00031 template<typename _Scalar, int NX=Dynamic, int NY=Dynamic>
00032 struct TestFunc1
00033 {
00034   typedef _Scalar Scalar;
00035   enum {
00036     InputsAtCompileTime = NX,
00037     ValuesAtCompileTime = NY
00038   };
00039   typedef Matrix<Scalar,InputsAtCompileTime,1> InputType;
00040   typedef Matrix<Scalar,ValuesAtCompileTime,1> ValueType;
00041   typedef Matrix<Scalar,ValuesAtCompileTime,InputsAtCompileTime> JacobianType;
00042 
00043   int m_inputs, m_values;
00044 
00045   TestFunc1() : m_inputs(InputsAtCompileTime), m_values(ValuesAtCompileTime) {}
00046   TestFunc1(int inputs, int values) : m_inputs(inputs), m_values(values) {}
00047 
00048   int inputs() const { return m_inputs; }
00049   int values() const { return m_values; }
00050 
00051   template<typename T>
00052   void operator() (const Matrix<T,InputsAtCompileTime,1>& x, Matrix<T,ValuesAtCompileTime,1>* _v) const
00053   {
00054     Matrix<T,ValuesAtCompileTime,1>& v = *_v;
00055 
00056     v[0] = 2 * x[0] * x[0] + x[0] * x[1];
00057     v[1] = 3 * x[1] * x[0] + 0.5 * x[1] * x[1];
00058     if(inputs()>2)
00059     {
00060       v[0] += 0.5 * x[2];
00061       v[1] += x[2];
00062     }
00063     if(values()>2)
00064     {
00065       v[2] = 3 * x[1] * x[0] * x[0];
00066     }
00067     if (inputs()>2 && values()>2)
00068       v[2] *= x[2];
00069   }
00070 
00071   void operator() (const InputType& x, ValueType* v, JacobianType* _j) const
00072   {
00073     (*this)(x, v);
00074 
00075     if(_j)
00076     {
00077       JacobianType& j = *_j;
00078 
00079       j(0,0) = 4 * x[0] + x[1];
00080       j(1,0) = 3 * x[1];
00081 
00082       j(0,1) = x[0];
00083       j(1,1) = 3 * x[0] + 2 * 0.5 * x[1];
00084 
00085       if (inputs()>2)
00086       {
00087         j(0,2) = 0.5;
00088         j(1,2) = 1;
00089       }
00090       if(values()>2)
00091       {
00092         j(2,0) = 3 * x[1] * 2 * x[0];
00093         j(2,1) = 3 * x[0] * x[0];
00094       }
00095       if (inputs()>2 && values()>2)
00096       {
00097         j(2,0) *= x[2];
00098         j(2,1) *= x[2];
00099 
00100         j(2,2) = 3 * x[1] * x[0] * x[0];
00101         j(2,2) = 3 * x[1] * x[0] * x[0];
00102       }
00103     }
00104   }
00105 };
00106 
00107 template<typename Func> void adolc_forward_jacobian(const Func& f)
00108 {
00109     typename Func::InputType x = Func::InputType::Random(f.inputs());
00110     typename Func::ValueType y(f.values()), yref(f.values());
00111     typename Func::JacobianType j(f.values(),f.inputs()), jref(f.values(),f.inputs());
00112 
00113     jref.setZero();
00114     yref.setZero();
00115     f(x,&yref,&jref);
00116 //     std::cerr << y.transpose() << "\n\n";;
00117 //     std::cerr << j << "\n\n";;
00118 
00119     j.setZero();
00120     y.setZero();
00121     AdolcForwardJacobian<Func> autoj(f);
00122     autoj(x, &y, &j);
00123 //     std::cerr << y.transpose() << "\n\n";;
00124 //     std::cerr << j << "\n\n";;
00125 
00126     VERIFY_IS_APPROX(y, yref);
00127     VERIFY_IS_APPROX(j, jref);
00128 }
00129 
00130 void test_forward_adolc()
00131 {
00132   adtl::ADOLC_numDir = NUMBER_DIRECTIONS;
00133 
00134   for(int i = 0; i < g_repeat; i++) {
00135     CALL_SUBTEST(( adolc_forward_jacobian(TestFunc1<double,2,2>()) ));
00136     CALL_SUBTEST(( adolc_forward_jacobian(TestFunc1<double,2,3>()) ));
00137     CALL_SUBTEST(( adolc_forward_jacobian(TestFunc1<double,3,2>()) ));
00138     CALL_SUBTEST(( adolc_forward_jacobian(TestFunc1<double,3,3>()) ));
00139     CALL_SUBTEST(( adolc_forward_jacobian(TestFunc1<double>(3,3)) ));
00140   }
00141 }


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