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
12template<typename _Scalar, int NX=Dynamic, int NY=Dynamic>
13struct 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
34struct 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
70void 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
91void 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
110void test_NumericalDiff()
111{
112    CALL_SUBTEST(test_forward());
113    CALL_SUBTEST(test_central());
114}
115