diff --git a/src/samplefunctional.hh b/src/samplefunctional.hh
index 5b8258b8efa4eee2a6f77cfebd6a8ab9709d6089..df53c3fdbb2e5e38325c062cf0388a4adfaea3ad 100644
--- a/src/samplefunctional.hh
+++ b/src/samplefunctional.hh
@@ -51,8 +51,12 @@ class SampleFunctional {
       return;
     } else if (pg * x >= 0 && mg * x >= 0) {
       ret = pg;
+      Dune::dverb << "## Expected directional derivative: " << -(pg * mg)
+                  << " (coordinates of the restriction)" << std::endl;
     } else if (pg * x <= 0 && mg * x <= 0) {
       ret = mg;
+      Dune::dverb << "## Expected directional derivative: " << -(mg * pg)
+                  << " (coordinates of the restriction)" << std::endl;
     } else {
       ret = project(smoothGradient(x), x);
     }
@@ -129,6 +133,15 @@ void minimise(const Functional J, const typename Functional::SmallVector x,
   MyDirectionalConvexFunctionType JRest(JRestA, JRestb, phi, x, descDir);
   // }}}
 
+  { // Debug
+    Interval<double> D;
+    JRest.subDiff(0, D);
+    Dune::dverb << "## Descent according to JRest.subDiff: " << D[1]
+                << std::endl;
+    assert(D[1] <=
+           0); // We should not be minimising in this direction otherwise
+  }
+
   // FIXME: default values are used
   Bisection bisection;
   int count;
diff --git a/src/test-gradient-method.cc b/src/test-gradient-method.cc
index 3e6c2d2c3a707a601cba3918f15e19648dcaa227..e12925add4eaa36674a10b4303ec77cd1bf21d68 100644
--- a/src/test-gradient-method.cc
+++ b/src/test-gradient-method.cc
@@ -108,6 +108,8 @@ void testSampleFunctionNonsmooth() {
     functionTester(J, start, 6);
   }
   std::cout << std::endl;
+  std::cout << std::endl;
+  std::cout << std::endl;
   {
     start = b;
     start /= (b.two_norm() - 1e-12); // Make sure the norm is above 1;