1 | // This file is part of Eigen, a lightweight C++ template library
|
---|
2 | // for linear algebra.
|
---|
3 | //
|
---|
4 | // Copyright (C) 2009 Thomas Capricelli <orzel@freehackers.org>
|
---|
5 |
|
---|
6 | #include <stdio.h>
|
---|
7 |
|
---|
8 | #include "main.h"
|
---|
9 | #include <unsupported/Eigen/NumericalDiff>
|
---|
10 |
|
---|
11 | // Generic functor
|
---|
12 | template<typename _Scalar, int NX=Dynamic, int NY=Dynamic>
|
---|
13 | struct Functor
|
---|
14 | {
|
---|
15 | typedef _Scalar Scalar;
|
---|
16 | enum {
|
---|
17 | InputsAtCompileTime = NX,
|
---|
18 | ValuesAtCompileTime = NY
|
---|
19 | };
|
---|
20 | typedef Matrix<Scalar,InputsAtCompileTime,1> InputType;
|
---|
21 | typedef Matrix<Scalar,ValuesAtCompileTime,1> ValueType;
|
---|
22 | typedef Matrix<Scalar,ValuesAtCompileTime,InputsAtCompileTime> JacobianType;
|
---|
23 |
|
---|
24 | int m_inputs, m_values;
|
---|
25 |
|
---|
26 | Functor() : m_inputs(InputsAtCompileTime), m_values(ValuesAtCompileTime) {}
|
---|
27 | Functor(int inputs, int values) : m_inputs(inputs), m_values(values) {}
|
---|
28 |
|
---|
29 | int inputs() const { return m_inputs; }
|
---|
30 | int values() const { return m_values; }
|
---|
31 |
|
---|
32 | };
|
---|
33 |
|
---|
34 | struct my_functor : Functor<double>
|
---|
35 | {
|
---|
36 | my_functor(void): Functor<double>(3,15) {}
|
---|
37 | int operator()(const VectorXd &x, VectorXd &fvec) const
|
---|
38 | {
|
---|
39 | double tmp1, tmp2, tmp3;
|
---|
40 | double y[15] = {1.4e-1, 1.8e-1, 2.2e-1, 2.5e-1, 2.9e-1, 3.2e-1, 3.5e-1,
|
---|
41 | 3.9e-1, 3.7e-1, 5.8e-1, 7.3e-1, 9.6e-1, 1.34, 2.1, 4.39};
|
---|
42 |
|
---|
43 | for (int i = 0; i < values(); i++)
|
---|
44 | {
|
---|
45 | tmp1 = i+1;
|
---|
46 | tmp2 = 16 - i - 1;
|
---|
47 | tmp3 = (i>=8)? tmp2 : tmp1;
|
---|
48 | fvec[i] = y[i] - (x[0] + tmp1/(x[1]*tmp2 + x[2]*tmp3));
|
---|
49 | }
|
---|
50 | return 0;
|
---|
51 | }
|
---|
52 |
|
---|
53 | int actual_df(const VectorXd &x, MatrixXd &fjac) const
|
---|
54 | {
|
---|
55 | double tmp1, tmp2, tmp3, tmp4;
|
---|
56 | for (int i = 0; i < values(); i++)
|
---|
57 | {
|
---|
58 | tmp1 = i+1;
|
---|
59 | tmp2 = 16 - i - 1;
|
---|
60 | tmp3 = (i>=8)? tmp2 : tmp1;
|
---|
61 | tmp4 = (x[1]*tmp2 + x[2]*tmp3); tmp4 = tmp4*tmp4;
|
---|
62 | fjac(i,0) = -1;
|
---|
63 | fjac(i,1) = tmp1*tmp2/tmp4;
|
---|
64 | fjac(i,2) = tmp1*tmp3/tmp4;
|
---|
65 | }
|
---|
66 | return 0;
|
---|
67 | }
|
---|
68 | };
|
---|
69 |
|
---|
70 | void test_forward()
|
---|
71 | {
|
---|
72 | VectorXd x(3);
|
---|
73 | MatrixXd jac(15,3);
|
---|
74 | MatrixXd actual_jac(15,3);
|
---|
75 | my_functor functor;
|
---|
76 |
|
---|
77 | x << 0.082, 1.13, 2.35;
|
---|
78 |
|
---|
79 | // real one
|
---|
80 | functor.actual_df(x, actual_jac);
|
---|
81 | // std::cout << actual_jac << std::endl << std::endl;
|
---|
82 |
|
---|
83 | // using NumericalDiff
|
---|
84 | NumericalDiff<my_functor> numDiff(functor);
|
---|
85 | numDiff.df(x, jac);
|
---|
86 | // std::cout << jac << std::endl;
|
---|
87 |
|
---|
88 | VERIFY_IS_APPROX(jac, actual_jac);
|
---|
89 | }
|
---|
90 |
|
---|
91 | void test_central()
|
---|
92 | {
|
---|
93 | VectorXd x(3);
|
---|
94 | MatrixXd jac(15,3);
|
---|
95 | MatrixXd actual_jac(15,3);
|
---|
96 | my_functor functor;
|
---|
97 |
|
---|
98 | x << 0.082, 1.13, 2.35;
|
---|
99 |
|
---|
100 | // real one
|
---|
101 | functor.actual_df(x, actual_jac);
|
---|
102 |
|
---|
103 | // using NumericalDiff
|
---|
104 | NumericalDiff<my_functor,Central> numDiff(functor);
|
---|
105 | numDiff.df(x, jac);
|
---|
106 |
|
---|
107 | VERIFY_IS_APPROX(jac, actual_jac);
|
---|
108 | }
|
---|
109 |
|
---|
110 | void test_NumericalDiff()
|
---|
111 | {
|
---|
112 | CALL_SUBTEST(test_forward());
|
---|
113 | CALL_SUBTEST(test_central());
|
---|
114 | }
|
---|