minor improvements to STIR Hessian code

- call accumulate_Hessian_times_input in the no-subset case, as this will now be a
bit faster (it's also a tiny bit less code)
- give correct name of parameters in cstir.h to avoid confusion
- add test for objective function with prior
- clarify RDP exclusion of test case

Author: KrisThielemans

src/xSTIR/cSTIR/include/sirf/STIR/cstir.h

@@ -173,9 +173,9 @@ extern "C" {
 	void* cSTIR_priorValue(void* ptr_p, void* ptr_i);
 	void* cSTIR_priorGradient(void* ptr_p, void* ptr_i);
     void* cSTIR_priorAccumulateHessianTimesInput
-        (void* ptr_prior, void* ptr_out, void* ptr_curr, void* ptr_inp);
+    (void* ptr_prior, void* ptr_curr, void* ptr_inp, void* ptr_out);
     void* cSTIR_priorComputeHessianTimesInput
-        (void* ptr_prior, void* ptr_out, void* ptr_cur, void* ptr_inp);
+        (void* ptr_prior, void* ptr_cur, void* ptr_inp, void* ptr_out);
 	void* cSTIR_computePriorGradient(void* ptr_p, void* ptr_i, void* ptr_g);
 	void* cSTIR_PLSPriorAnatomicalGradient(void* ptr_p, int dir);
 

src/xSTIR/cSTIR/include/sirf/STIR/stir_x.h

@@ -1188,11 +1188,8 @@ The actual algorithm is described in
 			output.fill(0.0);
 			if (subset >= 0)
 				accumulate_sub_Hessian_times_input(output, curr_image_est, input, subset);
-			else {
-				for (int s = 0; s < get_num_subsets(); s++) {
-					accumulate_sub_Hessian_times_input(output, curr_image_est, input, s);
-				}
-			}
+			else
+				accumulate_Hessian_times_input(output, curr_image_est, input);
 		}
 	};
 

src/xSTIR/pSTIR/tests/test_ObjectiveFunction.py

@@ -42,6 +42,9 @@ def setUp(self):
         am.set_up(templ,image)
         obj_fun = pet.make_Poisson_loglikelihood(acquired_data)
         obj_fun.set_acquisition_model(am)
+        prior = pet.QuadraticPrior()
+        prior.set_penalisation_factor(20)
+        obj_fun.set_prior(prior)
         obj_fun.set_up(image)
 
         self.obj_fun = obj_fun
@@ -62,14 +65,18 @@ def test_Hessian(self, subset=-1, eps=1e-3):
         """
         x = self.image
         dx = x.clone()
-        dx *= eps/dx.norm()
+        dx *= eps
         dx += eps/2
         y = x + dx
         gx = self.obj_fun.gradient(x, subset)
         gy = self.obj_fun.gradient(y, subset)
         dg = gy - gx
         Hdx = self.obj_fun.multiply_with_Hessian(x, dx, subset)
         q = (dg - Hdx).norm()/dg.norm()
+        print('norm of (x): %f' % x.norm())
+        print('norm of (x + dx): %f' % y.norm())
+        print('norm of grad(x): %f' % gx.norm())
+        print('norm of grad(x + dx): %f' % gy.norm())
         print('norm of grad(x + dx) - grad(x): %f' % dg.norm())
         print('norm of H(x)*dx: %f' % Hdx.norm())
         print('relative difference: %f' % q)

src/xSTIR/pSTIR/tests/tests_qp_lc_rdp.py

@@ -33,7 +33,7 @@ def Hessian_test(test, prior, x, eps=1e-3):
         """
     if x.norm() > 0:
         dx = x.clone()
-        dx *= eps/dx.norm()
+        dx *= eps
         dx += eps/2
     else:
         dx = x + eps
@@ -47,11 +47,17 @@ def Hessian_test(test, prior, x, eps=1e-3):
         q = (dg - Hdx).norm()/gynorm
     else:
         q = (dg - Hdx).norm()
-    #print('norm of grad(x + dx) - grad(x): %f' % dg.norm())
-    #print('norm of H(x)*dx: %f' % Hdx.norm())
-    #print('relative difference: %g' % q)
-    if dg.norm() == 0:
-        q = 0
+    # print('norm of x: %f, dx: %f' % (x.norm(), dx.norm()))
+    # print('norm of grad(x): %f, grad(x + dx): %f' % (gx.norm(), gy.norm()))
+    # print('norm of grad(x + dx) - grad(x): %f' % dg.norm())
+    # print('norm of H(x)*dx: %f' % Hdx.norm())
+    # print('relative difference: %g' % q)
+    if issubclass(type(prior), sirf.STIR.RelativeDifferencePrior):
+       if x.min() == x.max() and dx.min() == dx.max():
+          # skip test in this case, as grad(x+dx) = grad(x) = 0, but H dx is not (even analytically),
+          # although it is small.
+          # The difficult is knowing what the tolerance is for this test in that case
+          q = 0
     test.check_if_less(q, .01*eps)
 
 
@@ -67,6 +73,7 @@ def test_main(rec=False, verb=False, throw=True, no_ret_val=True):
     im_2 = im_0.get_uniform_copy(2)
 
     for im in [im_0, im_1, im_2]:
+      # print('-------------- new image (see test source) -----------------')
       for penalisation_factor in [0,1,4]:
         for kappa in [True, False]:
           priors = [sirf.STIR.QuadraticPrior(), sirf.STIR.LogcoshPrior(), sirf.STIR.RelativeDifferencePrior()]
@@ -97,6 +104,7 @@ def test_main(rec=False, verb=False, throw=True, no_ret_val=True):
 
             if isinstance(prior, sirf.STIR.RelativeDifferencePrior):
                 prior.set_epsilon(im.max()*.01)
+            # print(f"penalisation_factor {penalisation_factor}, kappa {kappa}, prior {type(prior).__name__}")
             Hessian_test(test, prior, im, 0.03)
 
     numpy.testing.assert_equal(test.failed, 0)