Home | History | Annotate | Download | only in test 
      1 // This file is part of Eigen, a lightweight C++ template library
      2 // for linear algebra.
      3 //
      4 // Copyright (C) 2009 Thomas Capricelli <orzel (at) 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 }
    115