add "increment" test for gradient (enabled for LM data)

also minor clean-up

src/include/stir/recon_buildblock/test/ObjectiveFunctionTests.h

@@ -15,7 +15,7 @@
   \brief Test skeleton for stir::GeneralisedObjectiveFunction and stir::GeneralisedPrior
 
   \author Kris Thielemans
-  \author Reobert Twyman Skelly
+  \author Robert Twyman Skelly
 */
 
 #include "stir/RunTests.h"
@@ -32,10 +32,6 @@ START_NAMESPACE_STIR
 
   This contains some numerical tests to check gradient and Hessian calculations.
 
-  Note that the gradient is computed numerically voxel by voxel which is obviously
-  terribly slow. A solution (for the test) would be to compute it only in
-  a subset of voxels or so. We'll leave this for later.
-
   Note that the test only works if the objective function is well-defined. For example,
   if certain projections are non-zero, while the model estimates them to be zero, the
   Poisson objective function is in theory infinite.
@@ -53,16 +49,22 @@ class ObjectiveFunctionTests : public RunTests
   typedef TargetT target_type;
 
   //! Test the gradient of the objective function by comparing to the numerical gradient via perturbation
-  /*! Note: \a target is non-\c const, as the code will add/subtract eps, but the actual values
+  /*!
+    If \a full_gradient=true, all elements in the gradient are tested (using single-element increments). This is slow.
+    Otherwise, the test checks that \f$ G dx \approx F(x+dx) - F(x) \f$.  dx is computed via construct_increment().
+
+    Note: \a target is non-\c const, as the code will add/subtract eps, but the actual values
     are not modified after the test exits.
   */
-  virtual Succeeded
-  test_gradient(const std::string& test_name, ObjectiveFunctionT& objective_function, TargetT& target, const float eps);
+  virtual Succeeded test_gradient(const std::string& test_name,
+                                  ObjectiveFunctionT& objective_function,
+                                  TargetT& target,
+                                  const float eps,
+                                  const bool full_gradient = true);
 
   //! Test the accumulate_Hessian_times_input of the objective function by comparing to the numerical result via perturbation
   /*!
-    This test checks that \f$ H dx \approx G(x+dx) - G(x) \f$.
-    \f$dx\f$ is currently computed by as <code>eps*(target / target.find_max() + 0.5)</code>.
+    This test checks that \f$ H dx \approx G(x+dx) - G(x) \f$. dx is computed via construct_increment().
   */
   virtual Succeeded
   test_Hessian(const std::string& test_name, ObjectiveFunctionT& objective_function, const TargetT& target, const float eps);
@@ -71,14 +73,23 @@ class ObjectiveFunctionTests : public RunTests
                                            ObjectiveFunctionT& objective_function,
                                            const TargetT& target,
                                            const float mult_factor = 1.F);
+
+protected:
+  //! Construct small increment for target
+  /*!
+    Result is <code>eps*(target / target.find_max() + 0.5)</code>, i.e. it is always
+    positive (if target is non-negative).
+  */
+  virtual shared_ptr<const TargetT> construct_increment(const TargetT& target, const float eps) const;
 };
 
 template <class ObjectiveFunctionT, class TargetT>
 Succeeded
 ObjectiveFunctionTests<ObjectiveFunctionT, TargetT>::test_gradient(const std::string& test_name,
                                                                    ObjectiveFunctionT& objective_function,
                                                                    TargetT& target,
-                                                                   const float eps)
+                                                                   const float eps,
+                                                                   const bool full_gradient)
 {
   shared_ptr<TargetT> gradient_sptr(target.get_empty_copy());
   shared_ptr<TargetT> gradient_2_sptr(target.get_empty_copy());
@@ -87,29 +98,47 @@ ObjectiveFunctionTests<ObjectiveFunctionT, TargetT>::test_gradient(const std::st
   this->set_tolerance(std::max(fabs(double(gradient_sptr->find_min())), double(gradient_sptr->find_max())) / 1000);
   info("Computing objective function at target");
   const double value_at_target = objective_function.compute_value(target);
-  auto target_iter = target.begin_all();
-  auto gradient_iter = gradient_sptr->begin_all();
-  auto gradient_2_iter = gradient_2_sptr->begin_all();
   bool testOK = true;
-  info("Computing gradient of objective function by numerical differences (this will take a while)");
-  while (target_iter != target.end_all())
+  if (full_gradient)
     {
-      *target_iter += eps;
-      const double value_at_inc = objective_function.compute_value(target);
-      *target_iter -= eps;
-      const float gradient_at_iter = static_cast<float>((value_at_inc - value_at_target) / eps);
-      *gradient_2_iter++ = gradient_at_iter;
-      testOK = testOK && this->check_if_equal(gradient_at_iter, *gradient_iter, "gradient");
-      ++target_iter;
-      ++gradient_iter;
+      info("Computing gradient of objective function by numerical differences (this will take a while)");
+      auto target_iter = target.begin_all();
+      auto gradient_iter = gradient_sptr->begin_all();
+      auto gradient_2_iter = gradient_2_sptr->begin_all();
+      while (target_iter != target.end_all())
+        {
+          *target_iter += eps;
+          const double value_at_inc = objective_function.compute_value(target);
+          *target_iter -= eps;
+          const float gradient_at_iter = static_cast<float>((value_at_inc - value_at_target) / eps);
+          *gradient_2_iter++ = gradient_at_iter;
+          testOK = testOK && this->check_if_equal(gradient_at_iter, *gradient_iter, "gradient");
+          ++target_iter;
+          ++gradient_iter;
+        }
+    }
+  else
+    {
+      /* test f(x+dx) - f(x) ~ dx^t G(x) */
+      shared_ptr<const TargetT> increment_sptr = this->construct_increment(target, eps);
+      shared_ptr<TargetT> target_plus_inc_sptr(target.clone());
+      *target_plus_inc_sptr += *increment_sptr;
+      const double value_at_inc = objective_function.compute_value(*target_plus_inc_sptr);
+      const double my_sum = std::inner_product(
+          gradient_sptr->begin_all_const(), gradient_sptr->end_all_const(), increment_sptr->begin_all_const(), 0.);
+
+      testOK = testOK && this->check_if_equal(value_at_inc - value_at_target, my_sum, "gradient");
     }
+
   if (!testOK)
     {
       std::cerr << "Numerical gradient test failed with for " + test_name + "\n";
-      std::cerr << "Writing diagnostic files " << test_name << "_target.hv, *gradient.hv, *numerical_gradient.hv\n";
+      std::cerr << "Writing diagnostic files " << test_name
+                << "_target.hv, *gradient.hv (and *numerical_gradient.hv if full gradient test is used)\n";
       write_to_file(test_name + "_target.hv", target);
       write_to_file(test_name + "_gradient.hv", *gradient_sptr);
-      write_to_file(test_name + "_numerical_gradient.hv", *gradient_2_sptr);
+      if (full_gradient)
+        write_to_file(test_name + "_numerical_gradient.hv", *gradient_2_sptr);
       return Succeeded::no;
     }
   else
@@ -118,20 +147,30 @@ ObjectiveFunctionTests<ObjectiveFunctionT, TargetT>::test_gradient(const std::st
     }
 }
 
+template <class ObjectiveFunctionT, class TargetT>
+shared_ptr<const TargetT>
+ObjectiveFunctionTests<ObjectiveFunctionT, TargetT>::construct_increment(const TargetT& target, const float eps) const
+{
+  shared_ptr<TargetT> increment_sptr(target.clone());
+  *increment_sptr *= eps / increment_sptr->find_max();
+  *increment_sptr += eps / 2;
+  return increment_sptr;
+}
+
 template <class ObjectiveFunctionT, class TargetT>
 Succeeded
 ObjectiveFunctionTests<ObjectiveFunctionT, TargetT>::test_Hessian(const std::string& test_name,
                                                                   ObjectiveFunctionT& objective_function,
                                                                   const TargetT& target,
                                                                   const float eps)
 {
+  info("Comparing Hessian*dx with difference of gradients");
+
   /* test G(x+dx) = G(x) + H dx + small stuff */
   shared_ptr<TargetT> gradient_sptr(target.get_empty_copy());
   shared_ptr<TargetT> gradient_2_sptr(target.get_empty_copy());
   shared_ptr<TargetT> output(target.get_empty_copy());
-  shared_ptr<TargetT> increment_sptr(target.clone());
-  *increment_sptr *= eps / increment_sptr->find_max();
-  *increment_sptr += eps / 2;
+  shared_ptr<const TargetT> increment_sptr = this->construct_increment(target, eps);
   shared_ptr<TargetT> target_plus_inc_sptr(target.clone());
   *target_plus_inc_sptr += *increment_sptr;
 
@@ -145,7 +184,6 @@ ObjectiveFunctionTests<ObjectiveFunctionT, TargetT>::test_Hessian(const std::str
   auto gradient_iter = gradient_sptr->begin_all_const();
   auto gradient_2_iter = gradient_2_sptr->begin_all_const();
   bool testOK = true;
-  info("Computing gradient of objective function by numerical differences (this will take a while)");
   while (output_iter != output->end_all())
     {
       testOK = testOK && this->check_if_equal(*gradient_2_iter - *gradient_iter, *output_iter, "Hessian*increment");
@@ -187,10 +225,9 @@ ObjectiveFunctionTests<ObjectiveFunctionT, TargetT>::test_Hessian_concavity(cons
   /// Compute dot(x,(H x))
   const float my_sum = std::inner_product(target.begin_all(), target.end_all(), output->begin_all(), 0.F) * mult_factor;
 
-  // test for a CONCAVE function
+  // test for a CONCAVE function (after multiplying with mult_factor of course)
   if (this->check_if_less(my_sum, 0))
     {
-      //    info("PASS: Computation of x^T H x = " + std::to_string(my_sum) + " < 0" (Hessian) and is therefore concave);
       return Succeeded::yes;
     }
   else
@@ -199,8 +236,8 @@ ObjectiveFunctionTests<ObjectiveFunctionT, TargetT>::test_Hessian_concavity(cons
       std::cerr << "FAIL: " + test_name + ": Computation of x^T H x = " + std::to_string(my_sum)
                 << " > 0 (Hessian) and is therefore NOT concave"
                 << "\n >target image max=" << target.find_max() << "\n >target image min=" << target.find_min()
-                << "\n >output image max=" << output->find_max() << "\n >output image min=" << output->find_min();
-      std::cerr << "Writing diagnostic files to " << test_name + "_concavity_out.hv etc\n.";
+                << "\n >output image max=" << output->find_max() << "\n >output image min=" << output->find_min() << '\n';
+      std::cerr << "Writing diagnostic files to " << test_name + "_concavity_out.hv, *target.hv\n";
       write_to_file(test_name + "_concavity_out.hv", *output);
       write_to_file(test_name + "_target.hv", target);
       return Succeeded::no;

src/recon_test/CMakeLists.txt

@@ -129,6 +129,7 @@ include(stir_test_exe_targets)
 # a test that uses MPI
 create_stir_mpi_test(test_PoissonLogLikelihoodWithLinearModelForMeanAndProjData.cxx "${STIR_LIBRARIES}" $<TARGET_OBJECTS:stir_registries>)
 
+# pass list-mode file as argument
 ADD_TEST(test_PoissonLogLikelihoodWithLinearModelForMeanAndListModeWithProjMatrixByBin
   test_PoissonLogLikelihoodWithLinearModelForMeanAndListModeWithProjMatrixByBin "${CMAKE_SOURCE_DIR}/recon_test_pack/PET_ACQ_small.l.hdr.STIR")
 

src/recon_test/test_PoissonLogLikelihoodWithLinearModelForMeanAndListModeWithProjMatrixByBin.cxx

@@ -119,15 +119,14 @@ void
 PoissonLogLikelihoodWithLinearModelForMeanAndListModeDataWithProjMatrixByBinTests::run_tests_for_objective_function(
     objective_function_type& objective_function, target_type& target)
 {
-  // TODO enable this, but far too slow ATM
-  // std::cerr << "----- testing Gradient\n";
-  // test_gradient("PoissonLLProjData", objective_function, target, 0.01F);
+  std::cerr << "----- testing Gradient\n";
+  test_gradient("PoissonLLListModeData", objective_function, target, 0.01F, /* full_gradient = */ false);
 
   std::cerr << "----- testing concavity via Hessian-vector product (accumulate_Hessian_times_input)\n";
-  test_Hessian_concavity("PoissonLLProjData", objective_function, target);
+  test_Hessian_concavity("PoissonLLListModeData", objective_function, target);
 
   std::cerr << "----- testing Hessian-vector product (accumulate_Hessian_times_input)\n";
-  test_Hessian("PoissonLLProjData", objective_function, target, 0.5F);
+  test_Hessian("PoissonLLListModeData", objective_function, target, 0.5F);
 }
 
 void

src/recon_test/test_PoissonLogLikelihoodWithLinearModelForMeanAndProjData.cxx

@@ -124,16 +124,16 @@ PoissonLogLikelihoodWithLinearModelForMeanAndProjDataTests::run_tests_for_object
     objective_function_type& objective_function, target_type& target)
 {
   std::cerr << "----- testing Gradient\n";
-  test_gradient("PoissonLLListModeData", objective_function, target, 0.01F);
+  test_gradient("PoissonLLProjData", objective_function, target, 0.01F);
 
   std::cerr << "----- testing concavity via Hessian-vector product (accumulate_Hessian_times_input)\n";
-  test_Hessian_concavity("PoissonLLListModeData", objective_function, target);
+  test_Hessian_concavity("PoissonLLProjData", objective_function, target);
 
   std::cerr << "----- testing approximate-Hessian-vector product (accumulate_Hessian_times_input)\n";
   test_approximate_Hessian_concavity(objective_function, target);
 
   std::cerr << "----- testing Hessian-vector product (accumulate_Hessian_times_input)\n";
-  test_Hessian("PoissonLLListModeData", objective_function, target, 0.5F);
+  test_Hessian("PoissonLLProjData", objective_function, target, 0.5F);
 
   if (!this->is_everything_ok())
     {