Implement testHessian method

dune/tnnmg/functionals/test/functionaltest.hh

@@ -11,8 +11,13 @@
 #include <vector>
 #include <cmath>
 #include <limits>
+#include <array>
+#include <iostream>
 
 #include <dune/common/exceptions.hh>
+#include <dune/common/classname.hh>
+
+#include <dune/istl/matrixindexset.hh>
 
 namespace Dune {
 
@@ -147,10 +152,82 @@ void testGradient(Functional& functional,
  * \param testPoints Test addHessian at these points
  */
 template <class Functional>
-void testHessian(const Functional& functional,
+void testHessian(Functional& functional,
                  const std::vector<typename Functional::VectorType>& testPoints)
 {
-  // TODO: Implement me!
+  for (auto&& testPoint : testPoints)
+  {
+    // Get the Hessian at the current test point as computed by 'functional'
+    MatrixIndexSet diagonalPattern;
+    diagonalPattern.resize(testPoint.size(), testPoint.size());
+    for (size_t i=0; i<testPoint.size(); i++)
+      diagonalPattern.add(i,i);
+
+    typename Functional::MatrixType hessian;
+    diagonalPattern.exportIdx(hessian);
+
+    hessian = 0;
+    functional.addHessian(testPoint, hessian);
+
+    // FD step size. Best value: fourth root of the machine precision
+    const double eps = std::pow(std::numeric_limits<double>::epsilon(),0.25);
+
+    for (size_t i=0; i<testPoint.size(); i++)
+    {
+      for (size_t j=0; j<testPoint.size(); j++)
+      {
+        for (size_t k=0; k<testPoint[i].size(); k++)
+        {
+          for (size_t l=0; l<testPoint[j].size(); l++)
+          {
+            // Do not compare anything if the functional claims to be non-differentiable
+            // at the test point in the given direction
+            functional.setVector(testPoint);
+            bool nonDifferentiableIK = functional.regularity(i, testPoint[i][k], k) > 1e10;
+            bool nonDifferentiableJL = functional.regularity(j, testPoint[j][l], l) > 1e10;
+            if (nonDifferentiableIK || nonDifferentiableJL)
+              continue;
+
+            // Compute an approximation to the derivative by finite differences
+            std::array<typename Functional::VectorType, 4> neighborPoints;
+            std::fill(neighborPoints.begin(), neighborPoints.end(), testPoint);
+
+            neighborPoints[0][i][k] += eps;
+            neighborPoints[0][j][l] += eps;
+            neighborPoints[1][i][k] -= eps;
+            neighborPoints[1][j][l] += eps;
+            neighborPoints[2][i][k] += eps;
+            neighborPoints[2][j][l] -= eps;
+            neighborPoints[3][i][k] -= eps;
+            neighborPoints[3][j][l] -= eps;
+
+            std::array<double,4> neighborValues;
+            for (int ii = 0; ii < 4; ii++)
+              neighborValues[ii] = functional(neighborPoints[ii]);
+
+            // The current partial derivative
+            double derivative = (i==j) ? hessian[i][j][k][l] : 0.0;
+
+            double finiteDiff = (neighborValues[0]
+                  - neighborValues[1] - neighborValues[2]
+                  + neighborValues[3]) / (4 * eps * eps);
+
+            // Check whether second derivative and its FD approximation match
+            if (std::abs(derivative - finiteDiff) > eps)
+            {
+              std::cout << "Bug in addHessian for functional type "
+                        << Dune::className<Functional>() << ":" << std::endl;
+              std::cout << "    Second derivative does not agree with FD approximation" << std::endl;
+              std::cout << "    Expected: " << derivative << ",  FD: " << finiteDiff << std::endl;
+              std::cout << std::endl;
+
+              DUNE_THROW(MathError,"");
+            }
+          }
+        }
+      }
+    }
+  }
 }
 
 /** \brief Test the directionalSubDiff method