Initial hessian computation

src/action.hpp

@@ -17,7 +17,7 @@
 
 template <typename Scalar, typename ScalarW = typename std::remove_const<Scalar>::type>
 inline Matrix<ScalarW,4,4>
-get_Q_matrix(array::ArrayView<Scalar, array::Shape<2,2,2> > Q)
+get_Q_matrix(const array::ArrayView<Scalar, array::Shape<2,2,2> >& Q)
 {
     auto g = Q.part(0u).matrix();
     auto gt = Q.part(1u).matrix();
@@ -40,7 +40,7 @@ get_Q_matrix(array::ArrayView<Scalar, array::Shape<2,2,2> > Q)
 }
 
 template <typename Scalar, typename Shape, typename Shape2, typename ScalarW = typename std::remove_const<Scalar>::type>
-inline ScalarW S_2(array::ArrayView<Scalar, Shape> Q, array::ArrayView<const Complex, Shape2> U, double h)
+inline ScalarW S_2(const array::ArrayView<Scalar, Shape>& Q, const array::ArrayView<const Complex, Shape2>& U, double h)
 {
     const int N_i[4] = {0, 1, -1, 0};
     const int N_j[4] = {1, 0, 0, -1};

src/core.cpp

@@ -41,7 +41,6 @@ public:
 
             Array<ADComplex, Shape<Dynamic,Dynamic,2,2,2>, std::vector<ADDouble> >
                 Q({nx_, ny_, 2, 2, 2});
-            Array<ADComplex, Shape<1>, std::vector<ADDouble> > res;
 
             auto a = Qval.reshape<Dynamic>({Qval.size()});
             auto b = Q.reshape<Dynamic>({Q.size()});
@@ -50,7 +49,20 @@ public:
 
             CppAD::Independent(Q.storage());
 
-            res(0u) = S_2(Q, U_, h_);
+            /*
+             * Take the real part only.
+             *
+             * Cauchy-Riemann guarantees the saddle point of real part is
+             * also saddle point of imaginary part.
+             *
+             * Non-analytic functions S(x,y) of N z=x+iy complex variables,
+             * do not generally have saddle points, because the gradient
+             * F = (d/dx Re S, d/dy Re S, d/dx Im S, d/dy Im S) contains
+             * more elements than there are variables.
+             */
+
+            Array<ADDouble, Shape<1> > res;
+            res(0u) = S_2(Q, U_, h_).real();
 
             return {Q.storage(), res.storage()};
         }
@@ -65,6 +77,61 @@ public:
             std::array<size_t, 1> shape = {jac.size()/2};
             return Array<Complex, Shape<Dynamic>, Vector<double> >(std::move(jac), shape);
         }
+
+    auto hess(QType Q)
+        {
+            Vector<double> x(reinterpret_cast<double *>(const_cast<Complex *>(Q.data())), 2 * Q.size());
+
+            CppAD::ADFun<double> f = eval_ad(Q);
+
+            /*
+             * Hessian of the real part
+             */
+
+            /*
+             * Compute sparsity pattern first
+             */
+            typedef Vector<size_t> SizeVector;
+            typedef Vector<double> DoubleVector;
+            typedef Vector<char> BoolVector;
+            typedef CppAD::sparse_rc<SizeVector> Sparsity;
+
+            /* Identity matrix */
+            size_t n      = 2*Q.size();
+            size_t nr     = n;
+            size_t nc     = n;
+            size_t nnz_in = n;
+            Sparsity pattern_in(nr, nc, nnz_in);
+            for(size_t k = 0; k < nnz_in; k++)
+            {
+                size_t r = k;
+                size_t c = k;
+                pattern_in.set(k, r, c);
+            }
+
+            /* Hessian sparsity: */
+            BoolVector s = {true};
+            Sparsity pattern;
+            f.rev_hes_sparsity(s, false, false, pattern);
+
+            /* Hessian */
+            CppAD::sparse_rcv<SizeVector, DoubleVector> subset(pattern);
+            CppAD::sparse_hes_work work;
+            DoubleVector w = {1.0};
+            std::string coloring = "cppad.symmetric";
+            size_t n_sweep = f.sparse_hes(x, w, subset, pattern, coloring, work);
+
+            if (n_sweep != 2)
+                throw std::runtime_error("hessian computation failed");
+
+            size_t nnz = subset.nnz();
+
+            Array<size_t, Shape<Dynamic>, SizeVector> row(std::move(SizeVector(subset.row())), {nnz});
+            Array<size_t, Shape<Dynamic>, SizeVector> col(std::move(SizeVector(subset.col())), {nnz});
+            Array<double, Shape<Dynamic>, DoubleVector> val(std::move(DoubleVector(subset.val())), {nnz});
+
+            return pybind11::make_tuple(row, col, val);
+        }
 };
 
 PYBIND11_MODULE(_core, m)