Implement minimise2() for 2D and test it

dune/tectonic/ellipticenergy.hh

@@ -52,6 +52,24 @@ template <int dim> class EllipticEnergy {
   }
 
   // returns false if the direction is tangential
+  void descentDirectionNew(SmallVector const &x, SmallVector &ret) const {
+    SmallVector pg;
+    upperGradient(x, pg);
+    SmallVector mg;
+    lowerGradient(x, mg);
+    double const pgx = pg * x;
+    double const mgx = mg * x;
+    if (pgx >= 0 && mgx >= 0) {
+      ret = pg;
+      ret *= -1;
+    } else if (pgx <= 0 && mgx <= 0) {
+      ret = mg;
+      ret *= -1;
+    } else {
+      assert(false);
+    }
+  }
+
   bool descentDirection(SmallVector const &x, SmallVector &ret) const {
     // Check the squared norm rather than each component because
     // complementaryProjection() divides by it
@@ -234,5 +252,65 @@ void minimise(Functional const &J, typename Functional::SmallVector &x,
     }
   }
 }
+
+// NOTE: only works for the 2D case
+template <class Functional>
+void minimise2(Functional const &J, typename Functional::SmallVector &x,
+               size_t steps, Bisection const &bisection) {
+  typedef typename Functional::SmallVector SmallVector;
+  minimisationInitialiser(J, bisection, x);
+
+  { // Make a single step where we already know we're not differentiable
+    SmallVector descDir;
+    if (x == SmallVector(0))
+      J.descentAtZero(descDir);
+    else
+      J.descentDirectionNew(x, descDir);
+    descentMinimisation(J, x, descDir, bisection);
+  }
+
+  // From here on, we're in a C1 region and stay there.
+  for (size_t i = 1; i < steps; ++i) {
+    SmallVector descDir;
+    J.descentDirectionNew(x, descDir);
+    descentMinimisation(J, x, descDir, bisection);
+  }
+}
+
+template <class Functional>
+void minimisationInitialiser(Functional const &J, Bisection const &bisection,
+                             typename Functional::SmallVector &startingPoint) {
+  typedef typename Functional::SmallVector SmallVector;
+
+  std::vector<double> const kinks = { 5, 10,
+                                      15 }; // FIXME: We happen to know these
+
+  SmallVector x_old(0);
+  double Jx_old = J(x_old);
+
+  for (auto &kink : kinks) {
+    SmallVector x1 = { kink, 0 }; // Random vector that has the right norm
+    SmallVector const descDir1 = { x1[1], -x1[0] };
+    tangentialMinimisation(J, x1, descDir1, bisection);
+    double const Jx1 = J(x1);
+
+    SmallVector x2 = { -x1[0], -x1[1] };
+    SmallVector const descDir2 = { x2[1], -x2[0] };
+    tangentialMinimisation(J, x2, descDir2, bisection);
+    double const Jx2 = J(x2);
+
+    double const Jx = (Jx2 < Jx1 ? Jx2 : Jx1);
+    SmallVector const x = (Jx2 < Jx1 ? x2 : x1);
+
+    if (Jx >= Jx_old) {
+      startingPoint = x_old;
+      return;
+    }
+
+    Jx_old = Jx;
+    x_old = x;
+  }
+  startingPoint = x_old;
+}
 }
 #endif

src/Makefile.am

@@ -13,7 +13,8 @@ check_PROGRAMS = \
 	test-gradient-sample-steep2 \
 	test-gradient-sample-verysteep \
 	test-gradient-sample2 \
-	test-gradient-trivial
+	test-gradient-trivial \
+	test-minimise2
 
 test_circle_1_SOURCES                      = test-circle.cc
 test_circle_1_CPPFLAGS                     = $(AM_CPPFLAGS) -DDUNE_TECTONIC_TEST_CIRCLE_SCALE=1
@@ -32,6 +33,7 @@ test_gradient_sample_steep2_SOURCES        = test-gradient-sample-steep2.cc
 test_gradient_sample_verysteep_SOURCES     = test-gradient-sample-verysteep.cc
 test_gradient_sample2_SOURCES              = test-gradient-sample2.cc
 test_gradient_trivial_SOURCES              = test-gradient-trivial.cc
+test_minimise2_SOURCES                     = test-minimise2.cc
 
 bin_PROGRAMS = \
 	one-body-sample-2D \

src/test-minimise2.cc

+/* Check that minimise2 always produces better results than minimise
+   when the latter starts from 0 */
+
+#ifdef HAVE_CONFIG_H
+#include "config.h"
+#endif
+
+#include <cassert>
+
+#include <boost/format.hpp>
+
+#include <dune/common/shared_ptr.hh>
+
+#include <dune/tectonic/ellipticenergy.hh>
+
+#include "test-gradient-method-nicefunction.hh"
+#include "test-gradient-method-helper.hh"
+
+int main() {
+  int const dim = 2;
+  typedef Dune::EllipticEnergy<dim> Functional;
+  typedef Functional::SmallMatrix SmallMatrix;
+  typedef Functional::SmallVector SmallVector;
+
+  SmallMatrix const A = { { 4, 1.5 }, { 1.5, 3 } };
+
+  std::vector<SmallVector> const bs = {
+    { 1, 2 }, { -8, -14 }, { -16, -28 }, { -24, -42 }, { -32, -56 },
+  };
+
+  auto const f = Dune::make_shared<Dune::ThreeKinkFunction const>();
+  auto const phi = Dune::make_shared<Functional::NonlinearityType const>(f);
+
+  std::vector<Functional> const Js = { Functional(A, bs[0], phi),
+                                       Functional(A, bs[1], phi),
+                                       Functional(A, bs[2], phi),
+                                       Functional(A, bs[3], phi),
+                                       Functional(A, bs[4], phi) };
+
+  Bisection const bisection(0.0, 1.0, 1e-12, false, 0);
+
+  size_t const runs = 5;
+
+  SmallVector x;
+
+  for (auto const &J : Js) {
+    x = { 0, 0 };
+    Dune::minimise(J, x, runs, bisection);
+    double const xmin = J(x);
+
+    x = { 0, 0 };
+    Dune::minimise2(J, x, runs, bisection);
+    double const xmin2 = J(x);
+    std::cout << xmin << std::endl;
+    std::cout << xmin2 << std::endl << std::endl;
+    assert(xmin2 <= xmin);
+  }
+}