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 | }
|
---|