try adept

meson.build

@@ -5,11 +5,15 @@ project('usadelndsoc', ['cpp'],
 
 cppad_dep = dependency('cppad')
 eigen3_dep = dependency('eigen3')
+adolc_dep = dependency('adolc')
 
 py_mod = import('python')
 py3 = py_mod.find_installation('python3')
 py3_dep = py3.dependency()
 
+cpp = meson.get_compiler('cpp')
+adept_dep = cpp.find_library('adept', has_headers : ['adept.h'])
+
 incdir_numpy = run_command(py3,
   [
     '-c',
@@ -38,7 +42,7 @@ pybind11_dep = declare_dependency(
   include_directories : [include_directories(incdir_pybind11)],
   dependencies : [py3_dep]
 )
-deps = [cppad_dep, eigen3_dep]
+deps = [cppad_dep, eigen3_dep, adept_dep, adolc_dep]
 executable('main', ['src/main.cpp'], dependencies : deps)
 
 subdir('usadelndsoc')

src/common.hpp

@@ -15,6 +15,8 @@
 #include <cppad/cppad.hpp>
 #include <cppad/example/cppad_eigen.hpp>
 
+#include <adept.h>
+
 #include "array.hpp"
 #include "vector.hpp"
 
@@ -22,18 +24,44 @@ typedef std::complex<double> Complex;
 typedef CppAD::AD<double> ADDouble;
 typedef std::complex<ADDouble> ADComplex;
 
-template <typename Scalar_, int Rows_, int Cols_>
-using Matrix = Eigen::Matrix<Scalar_, Rows_, Cols_, Eigen::RowMajor>;
-
-ADComplex adcomplex(const ADComplex&& cpx)
+class AdeptDouble
 {
-    return {cpx};
-}
+    typedef adept::adouble T;
+    
+public:
+    AdeptDouble() : v_() {}
+    AdeptDouble(const double& v) : v_(v) {}
+    AdeptDouble(const T& v) : v_(v) {}
+    AdeptDouble(const T&& v) : v_(std::move(v)) {}
 
-ADComplex adcomplex(const std::complex<double>&& cpx)
-{
-    return {cpx};
-}
+    T& value() { return v_; }
+    const T& value() const { return v_; }
+
+    AdeptDouble& operator+=(const AdeptDouble& other)
+        { v_ += other.v_; return *this; }
+    AdeptDouble& operator-=(const AdeptDouble& other)
+        { v_ -= other.v_; return *this; }
+    AdeptDouble& operator*=(const AdeptDouble& other)
+        { v_ *= other.v_; return *this; }
+    AdeptDouble& operator/=(const AdeptDouble& other)
+        { v_ /= other.v_; return *this; }
+
+    AdeptDouble operator-() const { return std::move(AdeptDouble(-v_)); }
+
+private:
+    adept::adouble v_;
+};
+
+AdeptDouble operator+(const AdeptDouble &a, const AdeptDouble &b) { return std::move(adept::adouble(a.value() + b.value())); }
+AdeptDouble operator-(const AdeptDouble &a, const AdeptDouble &b) { return std::move(adept::adouble(a.value() - b.value())); }
+AdeptDouble operator*(const AdeptDouble &a, const AdeptDouble &b) { return std::move(adept::adouble(a.value() * b.value())); }
+AdeptDouble operator/(const AdeptDouble &a, const AdeptDouble &b) { return std::move(adept::adouble(a.value() / b.value())); }
+
+
+typedef std::complex<AdeptDouble> AdeptComplex;
+
+template <typename Scalar_, int Rows_, int Cols_>
+using Matrix = Eigen::Matrix<Scalar_, Rows_, Cols_, Eigen::RowMajor>;
 
 namespace Eigen
 {
@@ -49,6 +77,16 @@ namespace Eigen
 }
 
 namespace CppAD {
+    ADComplex adcomplex(const ADComplex&& cpx)
+    {
+        return {cpx};
+    }
+
+    ADComplex adcomplex(const std::complex<double>&& cpx)
+    {
+        return {cpx};
+    }
+
     ADComplex operator*(const ADComplex& a, const double& b) { return a * adcomplex(b); }
     ADComplex operator*(const double& a, const ADComplex& b) { return adcomplex(a) * b; }
     ADComplex operator*(const ADComplex& a, const Complex& b) { return a * adcomplex(b); }
@@ -70,4 +108,47 @@ namespace CppAD {
     ADComplex operator/(const Complex& a, const ADComplex& b) { return adcomplex(a) / b; }
 }
 
+namespace Eigen
+{
+    /* Eigen binop maps: (double, ADDouble) -> ADDouble */
+    template<typename BinaryOp> struct ScalarBinaryOpTraits<AdeptDouble,double,BinaryOp> { typedef AdeptDouble ReturnType;  };
+    template<typename BinaryOp> struct ScalarBinaryOpTraits<double,AdeptDouble,BinaryOp> { typedef AdeptDouble ReturnType;  };
+
+    /* Eigen binop maps: (double/complex, AdeptComplex) -> AdeptComplex */
+    template<typename BinaryOp> struct ScalarBinaryOpTraits<AdeptComplex,double,BinaryOp> { typedef AdeptComplex ReturnType;  };
+    template<typename BinaryOp> struct ScalarBinaryOpTraits<double,AdeptComplex,BinaryOp> { typedef AdeptComplex ReturnType;  };
+    template<typename BinaryOp> struct ScalarBinaryOpTraits<AdeptComplex,Complex,BinaryOp> { typedef AdeptComplex ReturnType;  };
+    template<typename BinaryOp> struct ScalarBinaryOpTraits<Complex,AdeptComplex,BinaryOp> { typedef AdeptComplex ReturnType;  };
+}
+
+AdeptComplex adcomplex(const AdeptComplex&& cpx)
+{
+    return {cpx};
+}
+
+AdeptComplex adcomplex(const std::complex<double>&& cpx)
+{
+    return {cpx};
+}
+
+AdeptComplex operator*(const AdeptComplex& a, const double& b) { return a * adcomplex(b); }
+AdeptComplex operator*(const double& a, const AdeptComplex& b) { return adcomplex(a) * b; }
+AdeptComplex operator*(const AdeptComplex& a, const Complex& b) { return a * adcomplex(b); }
+AdeptComplex operator*(const Complex& a, const AdeptComplex& b) { return adcomplex(a) * b; }
+
+AdeptComplex operator+(const AdeptComplex& a, const double& b) { return a + adcomplex(b); }
+AdeptComplex operator+(const double& a, const AdeptComplex& b) { return adcomplex(a) + b; }
+AdeptComplex operator+(const AdeptComplex& a, const Complex& b) { return a + adcomplex(b); }
+AdeptComplex operator+(const Complex& a, const AdeptComplex& b) { return adcomplex(a) + b; }
+
+AdeptComplex operator-(const AdeptComplex& a, const double& b) { return a - adcomplex(b); }
+AdeptComplex operator-(const double& a, const AdeptComplex& b) { return adcomplex(a) - b; }
+AdeptComplex operator-(const AdeptComplex& a, const Complex& b) { return a - adcomplex(b); }
+AdeptComplex operator-(const Complex& a, const AdeptComplex& b) { return adcomplex(a) - b; }
+
+AdeptComplex operator/(const AdeptComplex& a, const double& b) { return a / adcomplex(b); }
+AdeptComplex operator/(const double& a, const AdeptComplex& b) { return adcomplex(a) / b; }
+AdeptComplex operator/(const AdeptComplex& a, const Complex& b) { return a / adcomplex(b); }
+AdeptComplex operator/(const Complex& a, const AdeptComplex& b) { return adcomplex(a) / b; }
+
 #endif

src/core.cpp

@@ -147,7 +147,7 @@ public:
             /* Hessian */
             DoubleVector w = {1.0};
             std::string coloring = "cppad.symmetric";
-            size_t n_sweep = f_.sparse_hes(x, w, hes_subset_, hes_pattern_, coloring, hes_work_);
+            f_.sparse_hes(x, w, hes_subset_, hes_pattern_, coloring, hes_work_);
 
             size_t nnz = hes_subset_.nnz();
 
@@ -159,6 +159,33 @@ public:
                 pybind11::make_tuple(std::move(val), pybind11::make_tuple(std::move(row), std::move(col))),
                 pybind11::make_tuple(x.size(), x.size()));
         }
+    auto grad_adept(const QType& Q)
+        {
+            if (Q.dim(0) != nx_ || Q.dim(1) != ny_)
+                throw std::out_of_range("Q array wrong size");
+
+            Vector<double> x = ravel(Q);
+            size_t n = x.size();
+
+            adept::Stack stack;
+
+            Array<AdeptComplex, Shape<Dynamic,Dynamic,2,2,2>, std::vector<AdeptDouble> >
+                Qad({nx_, ny_, 2, 2, 2});
+            Array<AdeptDouble, Shape<1> > res;
+
+            adept::set_values(reinterpret_cast<adept::adouble *>(Qad.storage().data()), n, x.data());
+
+            stack.new_recording();
+            res(0u) = S_2(Qad, U_, h_).real();
+
+            stack.independent(reinterpret_cast<adept::adouble *>(Qad.storage().data()), n);
+            stack.dependent(reinterpret_cast<adept::adouble *>(res.data()), 1);
+
+            Array<double, Shape<Dynamic> > jac({n});
+            stack.jacobian(jac.data());
+
+            return std::move(jac);
+        }
 };
 
 PYBIND11_MODULE(_core, m)
@@ -175,5 +202,6 @@ PYBIND11_MODULE(_core, m)
         .def("grad", &Action::grad)
         .def("hess", &Action::hess)
         .def("compute", &Action::compute)
+        .def("grad_adept", &Action::grad_adept)
         ;
 }