add a proof of concept autodiff jacobian helper class based on adolc

with unit test and FindAdolc cmake module

4 files changed, 202 insertions, 5 deletions

diff --git a/unsupported/CMakeLists.txt b/unsupported/CMakeLists.txt
index 6ebc94582..de3e0dc9f 100644
--- a/unsupported/CMakeLists.txt
+++ b/unsupported/CMakeLists.txt
@@ -2,7 +2,6 @@
 add_subdirectory(Eigen)
 
 if(EIGEN_BUILD_TESTS)
-  include(CTest)
   add_subdirectory(test)
 endif(EIGEN_BUILD_TESTS)
 
diff --git a/unsupported/Eigen/AdolcForward b/unsupported/Eigen/AdolcForward
index 2d65988d6..531b2ae5f 100644
--- a/unsupported/Eigen/AdolcForward
+++ b/unsupported/Eigen/AdolcForward
@@ -98,9 +98,9 @@ namespace adtl {
 
 }
 
-namespace Eigen { namespace unsupported { /*@}*/ } }
+namespace Eigen {
 
-template<> struct EigenNumTraits<adtl::adouble>
+template<> struct NumTraits<adtl::adouble>
 {
   typedef adtl::adouble Real;
   typedef adtl::adouble FloatingPoint;
@@ -113,4 +113,60 @@ template<> struct EigenNumTraits<adtl::adouble>
   };
 };
 
+template<typename Functor> class AdolcForwardJacobian : public Functor
+{
+  typedef adtl::adouble ActiveScalar;
+public:
+
+  AdolcForwardJacobian() : Functor() {}
+  AdolcForwardJacobian(const Functor& f) : Functor(f) {}
+
+  // forward constructors
+  template<typename T0>
+  AdolcForwardJacobian(const T0& a0) : Functor(a0) {}
+  template<typename T0, typename T1>
+  AdolcForwardJacobian(const T0& a0, const T1& a1) : Functor(a0, a1) {}
+  template<typename T0, typename T1, typename T2>
+  AdolcForwardJacobian(const T0& a0, const T1& a1, const T1& a2) : Functor(a0, a1, a2) {}
+
+  typedef typename Functor::InputType InputType;
+  typedef typename Functor::ValueType ValueType;
+  typedef typename Functor::JacobianType JacobianType;
+
+  typedef Matrix<ActiveScalar, InputType::SizeAtCompileTime, 1> ActiveInput;
+  typedef Matrix<ActiveScalar, ValueType::SizeAtCompileTime, 1> ActiveValue;
+
+  void operator() (const InputType& x, ValueType* v, JacobianType* _jac) const
+  {
+    ei_assert(v!=0);
+    if (!_jac)
+    {
+      Functor::operator()(x, v);
+      return;
+    }
+
+    JacobianType& jac = *_jac;
+
+    ActiveInput ax = x.template cast<ActiveScalar>();
+    ActiveValue av(jac.rows());
+
+    for (int j=0; j<jac.cols(); j++)
+      for (int i=0; i<jac.cols(); i++)
+        ax[i].setADValue(j, i==j ? 1 : 0);
+
+    Functor::operator()(ax, &av);
+
+    for (int i=0; i<jac.rows(); i++)
+    {
+      (*v)[i] = av[i].getValue();
+      for (int j=0; j<jac.cols(); j++)
+        jac.coeffRef(i,j) = av[i].getADValue(j);
+    }
+  }
+protected:
+
+};
+
+}
+
 #endif // EIGEN_ADLOC_FORWARD
diff --git a/unsupported/test/CMakeLists.txt b/unsupported/test/CMakeLists.txt
index 0efbd3d7d..c36bcad32 100644
--- a/unsupported/test/CMakeLists.txt
+++ b/unsupported/test/CMakeLists.txt
@@ -1,6 +1,16 @@
 
+include(EigenTesting)
+
 enable_testing()
 
-include(EigenTesting)
+find_package(Adolc)
+
+include_directories(../../test)
 
-# ei_add_test(foo "CXXFLAGS" "libs")
+if(ADOLC_FOUND)
+  include_directories(${ADOLC_INCLUDES})
+  ei_add_property(EIGEN_TESTED_BACKENDS "Adolc")
+  ei_add_test(forward_adolc " " ${ADOLC_LIBRARIES})
+else(ADOLC_FOUND)
+  ei_add_property(EIGEN_MISSING_BACKENDS "Adolc")
+endif(ADOLC_FOUND)
\ No newline at end of file
diff --git a/unsupported/test/forward_adolc.cpp b/unsupported/test/forward_adolc.cpp
new file mode 100644
index 000000000..016e20cdb
--- /dev/null
+++ b/unsupported/test/forward_adolc.cpp
@@ -0,0 +1,132 @@
+// This file is part of Eigen, a lightweight C++ template library
+// for linear algebra. Eigen itself is part of the KDE project.
+//
+// Copyright (C) 2008 Gael Guennebaud <g.gael@free.fr>
+//
+// Eigen is free software; you can redistribute it and/or
+// modify it under the terms of the GNU Lesser General Public
+// License as published by the Free Software Foundation; either
+// version 3 of the License, or (at your option) any later version.
+//
+// Alternatively, you can redistribute it and/or
+// modify it under the terms of the GNU General Public License as
+// published by the Free Software Foundation; either version 2 of
+// the License, or (at your option) any later version.
+//
+// Eigen is distributed in the hope that it will be useful, but WITHOUT ANY
+// WARRANTY; without even the implied warranty of MERCHANTABILITY or FITNESS
+// FOR A PARTICULAR PURPOSE. See the GNU Lesser General Public License or the
+// GNU General Public License for more details.
+//
+// You should have received a copy of the GNU Lesser General Public
+// License and a copy of the GNU General Public License along with
+// Eigen. If not, see <http://www.gnu.org/licenses/>.
+
+#include "main.h"
+#define NUMBER_DIRECTIONS 16
+#include <unsupported/Eigen/AdolcForward>
+
+int adtl::ADOLC_numDir;
+
+template<typename _Scalar, int NX, int NY>
+struct TestFunc1
+{
+  typedef _Scalar Scalar;
+  enum {
+    InputsAtCompileTime = NX,
+    ValuesAtCompileTime = NY
+  };
+  typedef Matrix<Scalar,InputsAtCompileTime,1> InputType;
+  typedef Matrix<Scalar,ValuesAtCompileTime,1> ValueType;
+  typedef Matrix<Scalar,ValuesAtCompileTime,InputsAtCompileTime> JacobianType;
+
+  template<typename T>
+  void operator() (const Matrix<T,InputsAtCompileTime,1>& x, Matrix<T,ValuesAtCompileTime,1>* _v) const
+  {
+    Matrix<T,ValuesAtCompileTime,1>& v = *_v;
+
+    v[0] = 2 * x[0] * x[0] + x[0] * x[1];
+    v[1] = 3 * x[1] * x[0] + 0.5 * x[1] * x[1];
+    if(NX>2)
+    {
+      v[0] += 0.5 * x[2];
+      v[1] += x[2];
+    }
+    if(NY>2)
+    {
+      v[2] = 3 * x[1] * x[0] * x[0];
+    }
+    if (NX>2 && NY>2)
+      v[2] *= x[2];
+  }
+
+  void operator() (const InputType& x, ValueType* v, JacobianType* _j) const
+  {
+    (*this)(x, v);
+
+    if(_j)
+    {
+      JacobianType& j = *_j;
+
+      j(0,0) = 4 * x[0] + x[1];
+      j(1,0) = 3 * x[1];
+
+      j(0,1) = x[0];
+      j(1,1) = 3 * x[0] + 2 * 0.5 * x[1];
+
+      if (NX>2)
+      {
+        j(0,2) = 0.5;
+        j(1,2) = 1;
+      }
+      if(NY>2)
+      {
+        j(2,0) = 3 * x[1] * 2 * x[0];
+        j(2,1) = 3 * x[0] * x[0];
+      }
+      if (NX>2 && NY>2)
+      {
+        j(2,0) *= x[2];
+        j(2,1) *= x[2];
+
+        j(2,2) = 3 * x[1] * x[0] * x[0];
+        j(2,2) = 3 * x[1] * x[0] * x[0];
+      }
+    }
+  }
+};
+
+template<typename Func> void adolc_forward_jacobian(const Func& f)
+{
+    typename Func::InputType x = Func::InputType::Random();
+    typename Func::ValueType y, yref;
+    typename Func::JacobianType j, jref;
+
+    jref.setZero();
+    yref.setZero();
+    f(x,&yref,&jref);
+//     std::cerr << y.transpose() << "\n\n";;
+//     std::cerr << j << "\n\n";;
+
+    j.setZero();
+    y.setZero();
+    AdolcForwardJacobian<Func> autoj(f);
+    autoj(x, &y, &j);
+//     std::cerr << y.transpose() << "\n\n";;
+//     std::cerr << j << "\n\n";;
+
+    VERIFY_IS_APPROX(y, yref);
+    VERIFY_IS_APPROX(j, jref);
+}
+
+void test_forward_adolc()
+{
+  adtl::ADOLC_numDir = NUMBER_DIRECTIONS;
+
+  for(int i = 0; i < g_repeat; i++) {
+    CALL_SUBTEST(( adolc_forward_jacobian(TestFunc1<double,2,2>()) ));
+    CALL_SUBTEST(( adolc_forward_jacobian(TestFunc1<double,2,3>()) ));
+    CALL_SUBTEST(( adolc_forward_jacobian(TestFunc1<double,3,2>()) ));
+    CALL_SUBTEST(( adolc_forward_jacobian(TestFunc1<double,3,3>()) ));
+  }
+}