Add unit test for Hessian via AutoDiffScalar

1 files changed, 50 insertions, 3 deletions

diff --git a/unsupported/test/autodiff.cpp b/unsupported/test/autodiff.cpp
index 087e7c542..1aa1b3d2d 100644
--- a/unsupported/test/autodiff.cpp
+++ b/unsupported/test/autodiff.cpp
@@ -129,6 +129,7 @@ template<typename Func> void forward_jacobian(const Func& f)
 
 
 // TODO also check actual derivatives!
+template <int>
 void test_autodiff_scalar()
 {
   Vector2f p = Vector2f::Random();
@@ -140,6 +141,7 @@ void test_autodiff_scalar()
 }
 
 // TODO also check actual derivatives!
+template <int>
 void test_autodiff_vector()
 {
   Vector2f p = Vector2f::Random();
@@ -153,6 +155,7 @@ void test_autodiff_vector()
   VERIFY_IS_APPROX(res.value(), foo(p));
 }
 
+template <int>
 void test_autodiff_jacobian()
 {
   CALL_SUBTEST(( forward_jacobian(TestFunc1<double,2,2>()) ));
@@ -162,12 +165,56 @@ void test_autodiff_jacobian()
   CALL_SUBTEST(( forward_jacobian(TestFunc1<double>(3,3)) ));
 }
 
+
+template <int>
+void test_autodiff_hessian()
+{
+  typedef AutoDiffScalar<VectorXd> AD;
+  typedef Matrix<AD,Eigen::Dynamic,1> VectorAD;
+  typedef AutoDiffScalar<VectorAD> ADD;
+  typedef Matrix<ADD,Eigen::Dynamic,1> VectorADD;
+  VectorADD x(2);
+  double s1 = internal::random<double>(), s2 = internal::random<double>(), s3 = internal::random<double>(), s4 = internal::random<double>();
+  x(0).value()=s1;
+  x(1).value()=s2;
+
+  //set unit vectors for the derivative directions (partial derivatives of the input vector)
+  x(0).derivatives().resize(2);
+  x(0).derivatives().setZero();
+  x(0).derivatives()(0)= 1;
+  x(1).derivatives().resize(2);
+  x(1).derivatives().setZero();
+  x(1).derivatives()(1)=1;
+
+  //repeat partial derivatives for the inner AutoDiffScalar
+  x(0).value().derivatives() = VectorXd::Unit(2,0);
+  x(1).value().derivatives() = VectorXd::Unit(2,1);
+
+  //set the hessian matrix to zero
+  for(int idx=0; idx<2; idx++) {
+      x(0).derivatives()(idx).derivatives()  = VectorXd::Zero(2);
+      x(1).derivatives()(idx).derivatives()  = VectorXd::Zero(2);
+  }
+
+  ADD y = sin(AD(s3)*x(0) + AD(s4)*x(1));
+
+  VERIFY_IS_APPROX(y.value().derivatives()(0), y.derivatives()(0).value());
+  VERIFY_IS_APPROX(y.value().derivatives()(1), y.derivatives()(1).value());
+  VERIFY_IS_APPROX(y.value().derivatives()(0), s3*std::cos(s1*s3+s2*s4));
+  VERIFY_IS_APPROX(y.value().derivatives()(1), s4*std::cos(s1*s3+s2*s4));
+  VERIFY_IS_APPROX(y.derivatives()(0).derivatives(), -std::sin(s1*s3+s2*s4)*Vector2d(s3*s3,s4*s3));
+  VERIFY_IS_APPROX(y.derivatives()(1).derivatives(),  -std::sin(s1*s3+s2*s4)*Vector2d(s3*s4,s4*s4));
+}
+
+
+
 void test_autodiff()
 {
   for(int i = 0; i < g_repeat; i++) {
-    CALL_SUBTEST_1( test_autodiff_scalar() );
-    CALL_SUBTEST_2( test_autodiff_vector() );
-    CALL_SUBTEST_3( test_autodiff_jacobian() );
+    CALL_SUBTEST_1( test_autodiff_scalar<1>() );
+    CALL_SUBTEST_2( test_autodiff_vector<1>() );
+    CALL_SUBTEST_3( test_autodiff_jacobian<1>() );
+    CALL_SUBTEST_4( test_autodiff_hessian<1>() );
   }
 }