Move minimisation code to a separate header

dune/tectonic/ellipticenergy.hh

@@ -6,11 +6,8 @@
 #include <dune/common/stdstreams.hh>
 
 #include <dune/fufem/interval.hh>
-#include <dune/tnnmg/problem-classes/bisection.hh>
 
-#include "mydirectionalconvexfunction.hh"
 #include "localnonlinearity.hh"
-#include "circularconvexfunction.hh"
 
 namespace Dune {
 template <int dim> class EllipticEnergy {
@@ -146,172 +143,5 @@ template <int dim> class EllipticEnergy {
     y.axpy(-dx / xx, x);
   }
 };
-
-template <class Functional>
-void descentMinimisation(Functional const &J,
-                         typename Functional::SmallVector &x,
-                         typename Functional::SmallVector const &descDir,
-                         Bisection const &bisection) {
-  typedef typename Functional::SmallVector SmallVector;
-  typedef typename Functional::NonlinearityType LocalNonlinearityType;
-
-  // {{{ Construct a restriction of J to the line x + t * descDir
-
-  /* We have
-
-     1/2 <A(u+xv),u+xv>-<b,u+xv> = 1/2 <Av,v> x^2 - <b-Au,v> x + <1/2 Au-b,u>
-
-     since A is symmetric.
-  */
-  SmallVector tmp = J.b;               //  b
-  J.A.mmv(x, tmp);                     //  b-Au
-  double const JRestb = tmp * descDir; // <b-Au,v>
-
-  J.A.mv(descDir, tmp);                //  Av
-  double const JRestA = tmp * descDir; // <Av,v>
-
-  MyDirectionalConvexFunction<LocalNonlinearityType> const JRest(
-      JRestA, JRestb, *J.phi, x, descDir);
-  // }}}
-
-  { // Debug
-    Interval<double> D;
-    JRest.subDiff(0, D);
-
-    dverb
-        << "## Directional derivative (as per subdifferential of restriction): "
-        << D[1] << " (coordinates of the restriction)" << std::endl;
-    /*
-      It is possible that this differs quite a lot from the
-      directional derivative computed in the descentDirection()
-      method:
-
-      If phi is nonsmooth at x, so that the directional
-      derivatives jump, and |x| is computed to be too small or too
-      large globally or locally, the locally computed
-      subdifferential and the globally computed subdifferential
-      will no longer coincide!
-
-      The assertion D[1] <= 0 may thus fail.
-    */
-  }
-
-  int count;
-  double const stepsize = bisection.minimize(JRest, 0.0, 0.0, count);
-  dverb << "Number of iterations in the bisection method: " << count
-        << std::endl;
-  ;
-
-  x.axpy(stepsize, descDir);
-}
-
-template <class Functional>
-void tangentialMinimisation(Functional const &J,
-                            typename Functional::SmallVector &x,
-                            typename Functional::SmallVector const &descDir,
-                            Bisection const &bisection) {
-  typedef typename Functional::SmallMatrix SmallMatrix;
-  typedef typename Functional::SmallVector SmallVector;
-
-  // We completely ignore the nonlinearity here -- when restricted
-  // to a circle, it just enters as a constant!
-  CircularConvexFunction<SmallMatrix, SmallVector> const JRest(J.A, J.b, x,
-                                                               descDir);
-
-  int count;
-  double const stepsize = bisection.minimize(JRest, 0.0, 1.0, count);
-  dverb << "Number of iterations in the bisection method: " << count
-        << std::endl;
-  ;
-
-  // Since x is used in the computation of the rhs, do not write to it directly
-  SmallVector tmp;
-  JRest.cartesian(stepsize, tmp);
-  x = tmp;
-}
-
-template <class Functional>
-void minimise(Functional const &J, typename Functional::SmallVector &x,
-              size_t steps, Bisection const &bisection) {
-  typedef typename Functional::SmallVector SmallVector;
-
-  for (size_t step = 0; step < steps; ++step) {
-    SmallVector descDir;
-    bool const linesearchp = J.descentDirection(x, descDir);
-
-    if (descDir == SmallVector(0.0))
-      return;
-
-    if (linesearchp) {
-      descentMinimisation(J, x, descDir, bisection);
-    } else {
-      Bisection slowBisection(bisection);
-      slowBisection.setFastQuadratic(false);
-
-      tangentialMinimisation(J, x, descDir, slowBisection);
-    }
-  }
-}
-
-// 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(0);
-    x2.axpy(-1, x1);
-    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

dune/tectonic/minimisation.hh

+#ifndef MINIMISATION_HH
+#define MINIMISATION_HH
+
+#include <dune/common/fmatrix.hh>
+#include <dune/common/fvector.hh>
+#include <dune/common/stdstreams.hh>
+
+#include <dune/fufem/interval.hh>
+#include <dune/tnnmg/problem-classes/bisection.hh>
+
+#include "circularconvexfunction.hh"
+#include "mydirectionalconvexfunction.hh"
+
+namespace Dune {
+template <class Functional>
+void descentMinimisation(Functional const &J,
+                         typename Functional::SmallVector &x,
+                         typename Functional::SmallVector const &descDir,
+                         Bisection const &bisection) {
+  typedef typename Functional::SmallVector SmallVector;
+  typedef typename Functional::NonlinearityType LocalNonlinearityType;
+
+  // {{{ Construct a restriction of J to the line x + t * descDir
+
+  /* We have
+
+     1/2 <A(u+xv),u+xv>-<b,u+xv> = 1/2 <Av,v> x^2 - <b-Au,v> x + <1/2 Au-b,u>
+
+     since A is symmetric.
+  */
+  SmallVector tmp = J.b;               //  b
+  J.A.mmv(x, tmp);                     //  b-Au
+  double const JRestb = tmp * descDir; // <b-Au,v>
+
+  J.A.mv(descDir, tmp);                //  Av
+  double const JRestA = tmp * descDir; // <Av,v>
+
+  MyDirectionalConvexFunction<LocalNonlinearityType> const JRest(
+      JRestA, JRestb, *J.phi, x, descDir);
+  // }}}
+
+  { // Debug
+    Interval<double> D;
+    JRest.subDiff(0, D);
+
+    dverb
+        << "## Directional derivative (as per subdifferential of restriction): "
+        << D[1] << " (coordinates of the restriction)" << std::endl;
+    /*
+      It is possible that this differs quite a lot from the
+      directional derivative computed in the descentDirection()
+      method:
+
+      If phi is nonsmooth at x, so that the directional
+      derivatives jump, and |x| is computed to be too small or too
+      large globally or locally, the locally computed
+      subdifferential and the globally computed subdifferential
+      will no longer coincide!
+
+      The assertion D[1] <= 0 may thus fail.
+    */
+  }
+
+  int count;
+  double const stepsize = bisection.minimize(JRest, 0.0, 0.0, count);
+  dverb << "Number of iterations in the bisection method: " << count
+        << std::endl;
+  ;
+
+  x.axpy(stepsize, descDir);
+}
+
+template <class Functional>
+void tangentialMinimisation(Functional const &J,
+                            typename Functional::SmallVector &x,
+                            typename Functional::SmallVector const &descDir,
+                            Bisection const &bisection) {
+  typedef typename Functional::SmallMatrix SmallMatrix;
+  typedef typename Functional::SmallVector SmallVector;
+
+  // We completely ignore the nonlinearity here -- when restricted
+  // to a circle, it just enters as a constant!
+  CircularConvexFunction<SmallMatrix, SmallVector> const JRest(J.A, J.b, x,
+                                                               descDir);
+
+  int count;
+  double const stepsize = bisection.minimize(JRest, 0.0, 1.0, count);
+  dverb << "Number of iterations in the bisection method: " << count
+        << std::endl;
+  ;
+
+  // Since x is used in the computation of the rhs, do not write to it directly
+  SmallVector tmp;
+  JRest.cartesian(stepsize, tmp);
+  x = tmp;
+}
+
+template <class Functional>
+void minimise(Functional const &J, typename Functional::SmallVector &x,
+              size_t steps, Bisection const &bisection) {
+  typedef typename Functional::SmallVector SmallVector;
+
+  for (size_t step = 0; step < steps; ++step) {
+    SmallVector descDir;
+    bool const linesearchp = J.descentDirection(x, descDir);
+
+    if (descDir == SmallVector(0.0))
+      return;
+
+    if (linesearchp) {
+      descentMinimisation(J, x, descDir, bisection);
+    } else {
+      Bisection slowBisection(bisection);
+      slowBisection.setFastQuadratic(false);
+
+      tangentialMinimisation(J, x, descDir, slowBisection);
+    }
+  }
+}
+
+// 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(0);
+    x2.axpy(-1, x1);
+    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

dune/tectonic/myblockproblem.hh

@@ -12,6 +12,7 @@
 
 #include "globalnonlinearity.hh"
 #include "localnonlinearity.hh"
+#include "minimisation.hh"
 #include "mydirectionalconvexfunction.hh"
 #include "ellipticenergy.hh"
 

src/test-gradient-method-helper.hh

@@ -6,6 +6,7 @@
 #include <dune/tnnmg/problem-classes/bisection.hh>
 
 #include <dune/tectonic/ellipticenergy.hh>
+#include <dune/tectonic/minimisation.hh>
 
 template <int dim>
 double functionTester(Dune::EllipticEnergy<dim> J,