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
FunctorFunctor26 Functor() : m_inputs(InputsAtCompileTime), m_values(ValuesAtCompileTime) {}
FunctorFunctor27 Functor(int inputs, int values) : m_inputs(inputs), m_values(values) {}
28
inputsFunctor29 int inputs() const { return m_inputs; }
valuesFunctor30 int values() const { return m_values; }
31
32 };
33
34 struct my_functor : Functor<double>
35 {
my_functormy_functor36 my_functor(void): Functor<double>(3,15) {}
operator ()my_functor37 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
actual_dfmy_functor53 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
test_forward()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
test_central()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
test_NumericalDiff()110 void test_NumericalDiff()
111 {
112 CALL_SUBTEST(test_forward());
113 CALL_SUBTEST(test_central());
114 }
115