[136] | 1 | // This file is part of Eigen, a lightweight C++ template library
|
---|
| 2 | // for linear algebra.
|
---|
| 3 | //
|
---|
| 4 | // Copyright (C) 2008 Gael Guennebaud <g.gael@free.fr>
|
---|
| 5 | //
|
---|
| 6 | // This Source Code Form is subject to the terms of the Mozilla
|
---|
| 7 | // Public License v. 2.0. If a copy of the MPL was not distributed
|
---|
| 8 | // with this file, You can obtain one at http://mozilla.org/MPL/2.0/.
|
---|
| 9 |
|
---|
| 10 | #include "main.h"
|
---|
| 11 | #include <Eigen/Dense>
|
---|
| 12 |
|
---|
| 13 | #define NUMBER_DIRECTIONS 16
|
---|
| 14 | #include <unsupported/Eigen/AdolcForward>
|
---|
| 15 |
|
---|
| 16 | int adtl::ADOLC_numDir;
|
---|
| 17 |
|
---|
| 18 | template<typename Vector>
|
---|
| 19 | EIGEN_DONT_INLINE typename Vector::Scalar foo(const Vector& p)
|
---|
| 20 | {
|
---|
| 21 | typedef typename Vector::Scalar Scalar;
|
---|
| 22 | return (p-Vector(Scalar(-1),Scalar(1.))).norm() + (p.array().sqrt().abs() * p.array().sin()).sum() + p.dot(p);
|
---|
| 23 | }
|
---|
| 24 |
|
---|
| 25 | template<typename _Scalar, int NX=Dynamic, int NY=Dynamic>
|
---|
| 26 | struct TestFunc1
|
---|
| 27 | {
|
---|
| 28 | typedef _Scalar Scalar;
|
---|
| 29 | enum {
|
---|
| 30 | InputsAtCompileTime = NX,
|
---|
| 31 | ValuesAtCompileTime = NY
|
---|
| 32 | };
|
---|
| 33 | typedef Matrix<Scalar,InputsAtCompileTime,1> InputType;
|
---|
| 34 | typedef Matrix<Scalar,ValuesAtCompileTime,1> ValueType;
|
---|
| 35 | typedef Matrix<Scalar,ValuesAtCompileTime,InputsAtCompileTime> JacobianType;
|
---|
| 36 |
|
---|
| 37 | int m_inputs, m_values;
|
---|
| 38 |
|
---|
| 39 | TestFunc1() : m_inputs(InputsAtCompileTime), m_values(ValuesAtCompileTime) {}
|
---|
| 40 | TestFunc1(int inputs, int values) : m_inputs(inputs), m_values(values) {}
|
---|
| 41 |
|
---|
| 42 | int inputs() const { return m_inputs; }
|
---|
| 43 | int values() const { return m_values; }
|
---|
| 44 |
|
---|
| 45 | template<typename T>
|
---|
| 46 | void operator() (const Matrix<T,InputsAtCompileTime,1>& x, Matrix<T,ValuesAtCompileTime,1>* _v) const
|
---|
| 47 | {
|
---|
| 48 | Matrix<T,ValuesAtCompileTime,1>& v = *_v;
|
---|
| 49 |
|
---|
| 50 | v[0] = 2 * x[0] * x[0] + x[0] * x[1];
|
---|
| 51 | v[1] = 3 * x[1] * x[0] + 0.5 * x[1] * x[1];
|
---|
| 52 | if(inputs()>2)
|
---|
| 53 | {
|
---|
| 54 | v[0] += 0.5 * x[2];
|
---|
| 55 | v[1] += x[2];
|
---|
| 56 | }
|
---|
| 57 | if(values()>2)
|
---|
| 58 | {
|
---|
| 59 | v[2] = 3 * x[1] * x[0] * x[0];
|
---|
| 60 | }
|
---|
| 61 | if (inputs()>2 && values()>2)
|
---|
| 62 | v[2] *= x[2];
|
---|
| 63 | }
|
---|
| 64 |
|
---|
| 65 | void operator() (const InputType& x, ValueType* v, JacobianType* _j) const
|
---|
| 66 | {
|
---|
| 67 | (*this)(x, v);
|
---|
| 68 |
|
---|
| 69 | if(_j)
|
---|
| 70 | {
|
---|
| 71 | JacobianType& j = *_j;
|
---|
| 72 |
|
---|
| 73 | j(0,0) = 4 * x[0] + x[1];
|
---|
| 74 | j(1,0) = 3 * x[1];
|
---|
| 75 |
|
---|
| 76 | j(0,1) = x[0];
|
---|
| 77 | j(1,1) = 3 * x[0] + 2 * 0.5 * x[1];
|
---|
| 78 |
|
---|
| 79 | if (inputs()>2)
|
---|
| 80 | {
|
---|
| 81 | j(0,2) = 0.5;
|
---|
| 82 | j(1,2) = 1;
|
---|
| 83 | }
|
---|
| 84 | if(values()>2)
|
---|
| 85 | {
|
---|
| 86 | j(2,0) = 3 * x[1] * 2 * x[0];
|
---|
| 87 | j(2,1) = 3 * x[0] * x[0];
|
---|
| 88 | }
|
---|
| 89 | if (inputs()>2 && values()>2)
|
---|
| 90 | {
|
---|
| 91 | j(2,0) *= x[2];
|
---|
| 92 | j(2,1) *= x[2];
|
---|
| 93 |
|
---|
| 94 | j(2,2) = 3 * x[1] * x[0] * x[0];
|
---|
| 95 | j(2,2) = 3 * x[1] * x[0] * x[0];
|
---|
| 96 | }
|
---|
| 97 | }
|
---|
| 98 | }
|
---|
| 99 | };
|
---|
| 100 |
|
---|
| 101 | template<typename Func> void adolc_forward_jacobian(const Func& f)
|
---|
| 102 | {
|
---|
| 103 | typename Func::InputType x = Func::InputType::Random(f.inputs());
|
---|
| 104 | typename Func::ValueType y(f.values()), yref(f.values());
|
---|
| 105 | typename Func::JacobianType j(f.values(),f.inputs()), jref(f.values(),f.inputs());
|
---|
| 106 |
|
---|
| 107 | jref.setZero();
|
---|
| 108 | yref.setZero();
|
---|
| 109 | f(x,&yref,&jref);
|
---|
| 110 | // std::cerr << y.transpose() << "\n\n";;
|
---|
| 111 | // std::cerr << j << "\n\n";;
|
---|
| 112 |
|
---|
| 113 | j.setZero();
|
---|
| 114 | y.setZero();
|
---|
| 115 | AdolcForwardJacobian<Func> autoj(f);
|
---|
| 116 | autoj(x, &y, &j);
|
---|
| 117 | // std::cerr << y.transpose() << "\n\n";;
|
---|
| 118 | // std::cerr << j << "\n\n";;
|
---|
| 119 |
|
---|
| 120 | VERIFY_IS_APPROX(y, yref);
|
---|
| 121 | VERIFY_IS_APPROX(j, jref);
|
---|
| 122 | }
|
---|
| 123 |
|
---|
| 124 | void test_forward_adolc()
|
---|
| 125 | {
|
---|
| 126 | adtl::ADOLC_numDir = NUMBER_DIRECTIONS;
|
---|
| 127 |
|
---|
| 128 | for(int i = 0; i < g_repeat; i++) {
|
---|
| 129 | CALL_SUBTEST(( adolc_forward_jacobian(TestFunc1<double,2,2>()) ));
|
---|
| 130 | CALL_SUBTEST(( adolc_forward_jacobian(TestFunc1<double,2,3>()) ));
|
---|
| 131 | CALL_SUBTEST(( adolc_forward_jacobian(TestFunc1<double,3,2>()) ));
|
---|
| 132 | CALL_SUBTEST(( adolc_forward_jacobian(TestFunc1<double,3,3>()) ));
|
---|
| 133 | CALL_SUBTEST(( adolc_forward_jacobian(TestFunc1<double>(3,3)) ));
|
---|
| 134 | }
|
---|
| 135 |
|
---|
| 136 | {
|
---|
| 137 | // simple instanciation tests
|
---|
| 138 | Matrix<adtl::adouble,2,1> x;
|
---|
| 139 | foo(x);
|
---|
| 140 | Matrix<adtl::adouble,Dynamic,Dynamic> A(4,4);;
|
---|
| 141 | A.selfadjointView<Lower>().eigenvalues();
|
---|
| 142 | }
|
---|
| 143 | }
|
---|