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