Use DirectionalConvexFunction<MyNonlinearityType>

where MyNonlinearityType is MyNonlinearity<dimension, Function>

In doing this,
 * Make directionalDerivative() assume |dir| = 1 and
 * Kill pseudoDirectionalDerivative

src/bisection-example-new.cc

@@ -38,52 +38,44 @@ class SampleFunctional {
     return y * v + func_(v.two_norm()); // <1/2 Av - b,v> + H(|v|)
   }
 
+  // Assumes |dir| = 1.
   double directionalDerivative(const SmallVector x,
                                const SmallVector dir) const {
-    double norm = dir.two_norm();
+    assert(abs(dir.two_norm() - 1) < 1e-8);
 
-    if (norm == 0)
-      DUNE_THROW(Dune::Exception, "Directional derivatives cannot be computed "
-                                  "w.r.t. the zero-direction.");
-
-    SmallVector tmp = dir;
-    tmp /= norm;
+    if (x == SmallVector(0.0))
+      return func_.rightDifferential(0);
 
-    return pseudoDirectionalDerivative(x, tmp);
+    return (x * dir > 0) ? PlusGrad(x) * dir : MinusGrad(x) * dir;
   }
 
   SmallVector minimise(const SmallVector x, unsigned int iterations) const {
     SmallVector descDir = ModifiedGradient(x);
 
+    if (descDir == SmallVector(0.0))
+      return SmallVector(0.0);
+
     /* 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 tmp2;
+    A_.mv(descDir, tmp2);
+    double const rest_A = tmp2 * descDir;
 
-    {
-      SmallVector tmp2;
-      A_.mv(descDir, tmp2);
-      double const rest_A = tmp2 * descDir;
+    SmallVector tmp3;
+    A_.mv(x, tmp3);
+    double const rest_b = (b_ - tmp3) * descDir;
 
-      SmallVector tmp3;
-      A_.mv(x, tmp3);
-      double const rest_b = (b_ - tmp3) * descDir;
+    typedef MyNonlinearity<dimension, Function> MyNonlinearityType;
+    MyNonlinearityType phi;
+    typedef DirectionalConvexFunction<MyNonlinearityType>
+    MyDirectionalConvexFunctionType;
+    MyDirectionalConvexFunctionType rest(rest_A, rest_b, phi, x, descDir);
 
-      typedef MyNonlinearity<dimension, SampleFunction> MyNonlinearityType;
-      MyNonlinearityType phi;
-      typedef DirectionalConvexFunction<MyNonlinearityType>
-      MyDirectionalConvexFunctionType;
-      MyDirectionalConvexFunctionType rest(rest_A, rest_b, phi, x, descDir);
-
-      // Experiment a bit
-      Interval<double> D;
-      rest.subDiff(0, D);
-    }
-
-    if (descDir == SmallVector(0.0))
-      return SmallVector(0.0);
+    Interval<double> D;
 
     // { Debugging
     Dune::dverb << "Starting at x with J(x) = " << operator()(x) << std::endl;
@@ -95,11 +87,9 @@ class SampleFunctional {
 
     double l = 0;
     double r = 1;
-    SmallVector tmp;
     while (true) {
-      tmp = x;
-      tmp.axpy(r, descDir);
-      if (pseudoDirectionalDerivative(tmp, descDir) >= 0)
+      rest.subDiff(r, D);
+      if (D[1] >= 0)
         break;
 
       l = r;
@@ -118,26 +108,26 @@ class SampleFunctional {
     }
 
     double m = l / 2 + r / 2;
-    SmallVector middle = SmallVector(0.0);
     for (unsigned int count = 0; count < iterations; ++count) {
-      Dune::dverb << "now at m = " << m << std::endl;
-      Dune::dverb << "Value of J here: " << operator()(x + middle) << std::endl;
-
-      middle = descDir;
-      middle *= m;
-
-      double pseudoDerivative =
-          pseudoDirectionalDerivative(x + middle, descDir);
-
-      if (pseudoDerivative < 0) {
+      { // Debugging
+        Dune::dverb << "now at m = " << m << std::endl;
+        SmallVector tmp = x;
+        tmp.axpy(m, descDir);
+        Dune::dverb << "Value of J here: " << operator()(tmp) << std::endl;
+      }
+
+      rest.subDiff(m, D);
+      if (D[1] < 0) {
         l = m;
         m = (m + r) / 2;
-      } else if (pseudoDerivative > 0) {
+      } else if (D[1] > 0) {
         r = m;
         m = (l + m) / 2;
       } else
         break;
     }
+    SmallVector middle = descDir;
+    middle *= m;
     return middle;
   }
 
@@ -167,20 +157,6 @@ class SampleFunctional {
     return y;
   }
 
-  // |dir|-times the directional derivative wrt dir/|dir|.  If only
-  // the sign of the directionalDerivative matters, this saves the
-  // cost of normalising.
-  double pseudoDirectionalDerivative(const SmallVector x,
-                                     const SmallVector dir) const {
-    if (x == SmallVector(0.0))
-      return func_.rightDifferential(0) * dir.two_norm();
-
-    if (x * dir > 0)
-      return PlusGrad(x) * dir;
-    else
-      return MinusGrad(x) * dir;
-  }
-
   SmallVector ModifiedGradient(const SmallVector x) const {
     if (x == SmallVector(0.0)) {
       SmallVector d = SmoothGrad(x);
@@ -234,14 +210,20 @@ void testSampleFunction() {
 
   SampleFunctional J(A, b);
 
-  std::cout << J.directionalDerivative(b, b) << std::endl;
-  assert(J.directionalDerivative(b, b) == 2 * sqrt(5) + 2);
+  SampleFunctional::SmallVector c = b;
+  c /= c.two_norm();
+
+  assert(J.directionalDerivative(b, c) == 2 * sqrt(5) + 2);
 
   SampleFunctional::SmallVector start = b;
   start *= 17;
   SampleFunctional::SmallVector correction = J.minimise(start, 20);
   assert(J(start + correction) <= J(start));
-  assert(std::abs(J(start + correction) + 0.254644) < 1e-8);
+
+  start += correction;
+  correction = J.minimise(start, 20);
+
+  assert(std::abs(J(start + correction) + 0.254631) < 1e-6);
   std::cout << "Arrived at J(...) = " << J(start + correction) << std::endl;
 }
 
@@ -260,8 +242,10 @@ void testTrivialFunction() {
 
   SampleFunctional J(A, b);
 
-  std::cout << J.directionalDerivative(b, b) << std::endl;
-  assert(J.directionalDerivative(b, b) == 2 * sqrt(5));
+  SampleFunctional::SmallVector c = b;
+  c /= c.two_norm();
+
+  assert(J.directionalDerivative(b, c) == 2 * sqrt(5));
 
   SampleFunctional::SmallVector start = b;
   start *= 17;