Fix an issue with the LL Hessian methods and adds tests

src/include/stir/recon_buildblock/PoissonLogLikelihoodWithLinearModelForMeanAndProjData.h

@@ -272,11 +272,11 @@ public  RegisteredParsingObject<PoissonLogLikelihoodWithLinearModelForMeanAndPro
     The Hessian (without penalty) is approximately given by:
     \f[ H_{jk} = - \sum_i P_{ij} h_i^{''}(y_i) P_{ik} \f]
     where
-    \f[ h_i(l) = y_i log (l) - l; h_i^{''}(y_i) = y_i / ((P \lambda)_i + a_i)^2; \f]
+    \f[ h_i(l) = y_i log (l) - l; h_i^{''}(y_i) = - y_i / ((P \lambda)_i + a_i)^2; \f]
     and \f$P_{ij} \f$ is the probability matrix.
 
     Hence
-    \f[ H_{jk} =  \sum_i P_{ij}(y_i / ((P \lambda)_i + a_i)^2) P_{ik} \f]
+    \f[ H_{jk} = - \sum_i P_{ij}(y_i / ((P \lambda)_i + a_i)^2) P_{ik} \f]
 
     This function is computationally expensive and can be approximated, see
     add_multiplication_with_approximate_sub_Hessian_without_penalty()

src/recon_buildblock/PoissonLogLikelihoodWithLinearModelForMeanAndProjData.cxx

@@ -1052,10 +1052,10 @@ actual_accumulate_sub_Hessian_times_input_without_penalty(TargetT& output,
 
   shared_ptr<TargetT> tmp(output.get_empty_copy());
   this->get_projector_pair().get_back_projector_sptr()->get_output(*tmp);
-  // output += tmp;
+  // output -= tmp;
   std::transform(output.begin_all(), output.end_all(),
                  tmp->begin_all(), output.begin_all(),
-                 std::plus<typename TargetT::full_value_type>());
+                 std::minus<typename TargetT::full_value_type>());
 
   return Succeeded::yes;
 }

src/recon_test/test_PoissonLogLikelihoodWithLinearModelForMeanAndProjData.cxx

@@ -107,6 +107,14 @@ class PoissonLogLikelihoodWithLinearModelForMeanAndProjDataTests : public RunTes
   /*! Note that this function is not specific to PoissonLogLikelihoodWithLinearModelForMeanAndProjData */
   void run_tests_for_objective_function(GeneralisedObjectiveFunction<target_type>& objective_function,
                                         target_type& target);
+
+  //! Test the gradient of the objective function by comparing to the numerical gradient via perturbation
+  void test_objective_function_gradient(GeneralisedObjectiveFunction<target_type>& objective_function,
+                                        target_type& target);
+
+  //! Test the Hessian of the objective function by testing the (x^T Hx > 0) condition
+  void test_objective_function_Hessian_concavity(GeneralisedObjectiveFunction<target_type>& objective_function,
+                                                 target_type& target);
 };
 
 PoissonLogLikelihoodWithLinearModelForMeanAndProjDataTests::
@@ -118,6 +126,20 @@ void
 PoissonLogLikelihoodWithLinearModelForMeanAndProjDataTests::
 run_tests_for_objective_function(GeneralisedObjectiveFunction<PoissonLogLikelihoodWithLinearModelForMeanAndProjDataTests::target_type>& objective_function,
                                  PoissonLogLikelihoodWithLinearModelForMeanAndProjDataTests::target_type& target) {
+  std::cerr << "----- testing Gradient\n";
+  test_objective_function_gradient(objective_function,
+                                   target);
+
+  std::cerr << "----- testing Hessian-vector product (accumulate_Hessian_times_input)\n";
+  test_objective_function_Hessian_concavity(objective_function,
+                                            target);
+}
+
+void
+PoissonLogLikelihoodWithLinearModelForMeanAndProjDataTests::
+test_objective_function_gradient(GeneralisedObjectiveFunction<target_type> &objective_function,
+                                 target_type &target)
+{
   shared_ptr<target_type> gradient_sptr(target.get_empty_copy());
   shared_ptr<target_type> gradient_2_sptr(target.get_empty_copy());
   const int subset_num = 0;
@@ -166,6 +188,38 @@ run_tests_for_objective_function(GeneralisedObjectiveFunction<PoissonLogLikeliho
     }
 
 }
+void
+PoissonLogLikelihoodWithLinearModelForMeanAndProjDataTests::
+test_objective_function_Hessian_concavity(GeneralisedObjectiveFunction<target_type> &objective_function,
+                                          target_type &target){
+  /// setup images
+  shared_ptr<target_type> output(target.get_empty_copy());
+
+  /// Compute H x
+  objective_function.accumulate_Hessian_times_input(*output, target, target);
+
+  /// Compute dot(x,(H x))
+  float my_sum = 0.0;
+  {
+    target_type::full_iterator target_iter = target.begin_all();
+    target_type::full_iterator output_iter = output->begin_all();
+    while (target_iter != target.end_all())// && testOK)
+    {
+      my_sum += *target_iter * *output_iter;
+      ++target_iter; ++output_iter;
+    }
+  }
+  // test for a CONCAVE function
+  if (this->check_if_less( my_sum, 0)) {
+//    info("PASS: Computation of x^T H x = " + std::to_string(my_sum) + " < 0" and is therefore concave);
+  } else {
+    // print to console the FAILED configuration
+    info("FAIL: Computation of x^T H x = " + std::to_string(my_sum) + " > 0 and is therefore NOT concave" +
+         "\n >target image max=" + std::to_string(target.find_max()) +
+         "\n >target image min=" + std::to_string(target.find_min()));
+  }
+}
+
 
 void
 PoissonLogLikelihoodWithLinearModelForMeanAndProjDataTests::