Use bisection class from dune-tnnmg

src/bisection-example-new.cc

@@ -49,7 +49,7 @@ class SampleFunctional {
     return (x * dir > 0) ? PlusGrad(x) * dir : MinusGrad(x) * dir;
   }
 
-  SmallVector minimise(const SmallVector x, unsigned int iterations) const {
+  SmallVector minimise(const SmallVector x) const {
     SmallVector descDir = ModifiedGradient(x);
 
     if (descDir == SmallVector(0.0))
@@ -75,57 +75,14 @@ class SampleFunctional {
     MyDirectionalConvexFunctionType;
     MyDirectionalConvexFunctionType rest(rest_A, rest_b, phi, x, descDir);
 
-    Interval<double> D;
-
-    // { Debugging
-    Dune::dverb << "Starting at x with J(x) = " << operator()(x) << std::endl;
-    SmallVector debug_normalised_dir = descDir;
-    debug_normalised_dir /= debug_normalised_dir.two_norm();
-    Dune::dverb << "Minimizing in direction w with dJ(x,w) = "
-                << directionalDerivative(x, debug_normalised_dir) << std::endl;
-    // }
-
-    double l = 0;
-    double r = 1;
-    while (true) {
-      rest.subDiff(r, D);
-      if (D[1] >= 0)
-        break;
-
-      l = r;
-      r *= 2;
-      Dune::dverb << "Widened interval!" << std::endl;
-    }
-    Dune::dverb << "Interval now [" << l << "," << r << "]" << std::endl;
-
-    { // Debugging
-      SmallVector tmpl = x;
-      tmpl.axpy(l, descDir);
-      SmallVector tmpr = x;
-      tmpr.axpy(r, descDir);
-      assert(directionalDerivative(tmpl, debug_normalised_dir) < 0);
-      assert(directionalDerivative(tmpr, debug_normalised_dir) > 0);
-    }
+    Bisection bisection; // FIXME: default values
+    int count;
+    // FIXME: does x_old = 1 make any sense?!
+    double const m = bisection.minimize(rest, 0.0, 1.0, count);
+    Dune::dverb << "Number of iterations in the bisection method: " << count
+                << std::endl;
+    ;
 
-    double m = l / 2 + r / 2;
-    for (unsigned int count = 0; count < iterations; ++count) {
-      { // 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 (D[1] > 0) {
-        r = m;
-        m = (l + m) / 2;
-      } else
-        break;
-    }
     SmallVector middle = descDir;
     middle *= m;
     return middle;
@@ -220,8 +177,9 @@ void testSampleFunction() {
   SampleFunctional::SmallVector correction;
 
   for (int i = 1; i <= 6; ++i) {
-    correction = J.minimise(start, 20);
+    correction = J.minimise(start);
     start += correction;
+    std::cout << "New value: J(...) = " << J(start) << std::endl;
   }
   std::cout << "Arrived at J(...) = " << J(start) << std::endl;
 }
@@ -251,8 +209,9 @@ void testTrivialFunction() {
   SampleFunctional::SmallVector correction;
 
   for (int i = 1; i <= 5; ++i) {
-    correction = J.minimise(start, 20);
+    correction = J.minimise(start);
     start += correction;
+    std::cout << "New value: J(...) = " << J(start) << std::endl;
   }
   std::cout << "Arrived at J(...) = " << J(start) << std::endl;
 }