Computing Hessian products with image data objects and related stuff

examples/Python/PET/osem_reconstruction.py

@@ -128,6 +128,8 @@ def main():
     print('setting up, please wait...')
     recon.set_up(image)
 
+    obj_fun.test_Hessian(image, subset=0)
+
     # set the initial image estimate
     recon.set_current_estimate(image)
 

src/xSTIR/cSTIR/cstir.cpp

@@ -1,6 +1,6 @@
 /*
 SyneRBI Synergistic Image Reconstruction Framework (SIRF)
-Copyright 2017 - 2023 Rutherford Appleton Laboratory STFC
+Copyright 2017 - 2024 Rutherford Appleton Laboratory STFC
 Copyright 2018 - 2024 University College London.
 
 This is software developed for the Collaborative Computational
@@ -417,12 +417,28 @@ void* cSTIR_objectFromFile(const char* name, const char* filename)
 	CATCH;
 }
 
+typedef xSTIR_PoissonLLhLinModMeanListDataProjMatBin3DF LMObjFun;
+
+extern "C"
+void* cSTIR_objFunListModeSetInterval(void* ptr_f, size_t ptr_data)
+{
+	try {
+		auto objFun = objectFromHandle<LMObjFun>(ptr_f);
+		// Does not compile!!!
+	    //xSTIR_PoissonLLhLinModMeanListDataProjMatBin3DF& objFun =
+		//	objectFromHandle<xSTIR_PoissonLLhLinModMeanListDataProjMatBin3D>(ptr_f);
+		float* data = (float*)ptr_data;
+		objFun.set_time_interval((double)data[0], (double)data[1]);
+		return (void*)new DataHandle;
+	}
+	CATCH;
+}
+
 extern "C"
 void* cSTIR_setListmodeToSinogramsInterval(void* ptr_lm2s, size_t ptr_data)
 {
 	try {
-		ListmodeToSinograms& lm2s = 
-			objectFromHandle<ListmodeToSinograms>(ptr_lm2s);
+		auto lm2s = objectFromHandle<ListmodeToSinograms>(ptr_lm2s);
 		float *data = (float *)ptr_data;
 		lm2s.set_time_interval((double)data[0], (double)data[1]);
 		return (void*)new DataHandle;
@@ -1249,6 +1265,51 @@ cSTIR_computeObjectiveFunctionGradientNotDivided(void* ptr_f, void* ptr_i, int s
 	CATCH;
 }
 
+extern "C"
+void*
+cSTIR_objectiveFunctionAccumulateHessianTimesInput
+    (void* ptr_fun, void* ptr_est, void* ptr_inp, int subset, void* ptr_out)
+{
+	try {
+		auto& fun = objectFromHandle<ObjectiveFunction3DF>(ptr_fun);
+		auto& est = objectFromHandle<STIRImageData>(ptr_est);
+		auto& inp = objectFromHandle<STIRImageData>(ptr_inp);
+		auto& out = objectFromHandle<STIRImageData>(ptr_out);
+		auto& curr_est = est.data();
+		auto& input    = inp.data();
+		auto& output   = out.data();
+		output.fill(0.0);
+		if (subset >= 0)
+			fun.accumulate_sub_Hessian_times_input(output, curr_est, input, subset);
+		else {
+			for (int s = 0; s < fun.get_num_subsets(); s++) {
+				fun.accumulate_sub_Hessian_times_input(output, curr_est, input, s);
+			}
+		}
+		return (void*) new DataHandle;
+	}
+	CATCH;
+}
+
+extern "C"
+void*
+cSTIR_objectiveFunctionComputeHessianTimesInput
+    (void* ptr_fun, void* ptr_est, void* ptr_inp, int subset, void* ptr_out)
+{
+	try {
+		auto& fun = objectFromHandle<xSTIR_GeneralisedObjectiveFunction3DF>(ptr_fun);
+		auto& est = objectFromHandle<STIRImageData>(ptr_est);
+		auto& inp = objectFromHandle<STIRImageData>(ptr_inp);
+		auto& out = objectFromHandle<STIRImageData>(ptr_out);
+		auto& curr_est = est.data();
+		auto& input    = inp.data();
+		auto& output   = out.data();
+		fun.multiply_with_Hessian(output, curr_est, input, subset);
+		return (void*) new DataHandle;
+	}
+	CATCH;
+}
+
 extern "C"
 void*
 cSTIR_setupPrior(void* ptr_p, void* ptr_i)
@@ -1288,7 +1349,7 @@ void*
 cSTIR_priorGradient(void* ptr_p, void* ptr_i)
 {
 	try {
-		Prior3DF& prior = objectFromHandle<Prior3DF>(ptr_p);
+		Prior3DF& prior = objectFromHandle<stir::GeneralisedPrior <Image3DF> >(ptr_p);
 		STIRImageData& id = objectFromHandle<STIRImageData>(ptr_i);
 		Image3DF& image = id.data();
 		shared_ptr<STIRImageData> sptr(new STIRImageData(image));
@@ -1299,6 +1360,24 @@ cSTIR_priorGradient(void* ptr_p, void* ptr_i)
 	CATCH;
 }
 
+extern "C"
+void*
+cSTIR_priorAccumulateHessianTimesInput(void* ptr_prior, void* ptr_out, void* ptr_cur, void* ptr_inp)
+{
+	try {
+		auto& prior = objectFromHandle<stir::GeneralisedPrior <Image3DF> >(ptr_prior);
+		auto& out = objectFromHandle<STIRImageData>(ptr_out);
+		auto& cur = objectFromHandle<STIRImageData>(ptr_cur);
+		auto& inp = objectFromHandle<STIRImageData>(ptr_inp);
+		auto& output  = out.data();
+		auto& current = cur.data();
+		auto& input   = inp.data();
+		prior.accumulate_Hessian_times_input(output, current, input);
+		return (void*) new DataHandle;
+	}
+	CATCH;
+}
+
 extern "C"
 void*
 cSTIR_computePriorGradient(void* ptr_p, void* ptr_i, void* ptr_g)

src/xSTIR/cSTIR/cstir_p.cpp

@@ -1012,6 +1012,11 @@ sirf::cSTIR_setPoissonLogLikelihoodWithLinearModelForMeanAndListModeDataWithProj
     else if (sirf::iequals(name, "cache_max_size")) {
         obj_fun.set_cache_max_size(dataFromHandle<int>(hv));
     }
+    else if (sirf::iequals(name, "subsensitivity_filenames"))
+    {
+        std::string s(charDataFromDataHandle(hv));
+        obj_fun.set_subsensitivity_filenames(s.c_str());
+    }
     else
         return parameterNotFound(name, __FILE__, __LINE__);
     return new DataHandle;

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

@@ -65,6 +65,7 @@ extern "C" {
 	void* cSTIR_convertListmodeToSinograms(void* ptr);
 	void* cSTIR_computeRandoms(void* ptr);
     void* cSTIR_lm_num_prompts_exceeds_threshold(void* ptr, const float threshold);
+    void* cSTIR_objFunListModeSetInterval(void* ptr_f, size_t ptr_data);
 
 	// Data processor methods
 	void* cSTIR_setupImageDataProcessor(const void* ptr_p, void* ptr_i);
@@ -153,11 +154,17 @@ extern "C" {
 		(void* ptr_f, void* ptr_i, int subset);
     void* cSTIR_computeObjectiveFunctionGradientNotDivided
         (void* ptr_f, void* ptr_i, int subset, void* ptr_g);
+    void* cSTIR_objectiveFunctionAccumulateHessianTimesInput
+        (void* ptr_fun, void* ptr_est, void* ptr_inp, int subset, void* ptr_out);
+    void* cSTIR_objectiveFunctionComputeHessianTimesInput
+        (void* ptr_fun, void* ptr_est, void* ptr_inp, int subset, void* ptr_out);
 
 	// Prior methods
 	void* cSTIR_setupPrior(void* ptr_p, void* ptr_i);
 	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* 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

@@ -1067,6 +1067,19 @@ The actual algorithm is described in
 	class xSTIR_GeneralisedObjectiveFunction3DF :
 		public stir::GeneralisedObjectiveFunction < Image3DF > {
 	public:
+		void multiply_with_Hessian(Image3DF& output, const Image3DF& curr_image_est,
+            const Image3DF& input, const int subset) const
+		{
+			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);
+				}
+			}
+		}
+
 //		bool post_process() {
 //			return post_processing();
 //		}
@@ -1118,6 +1131,15 @@ The actual algorithm is described in
                  set_cache_path(filepath);
         }
 
+        void set_time_interval(double start, double stop)
+        {
+	        std::pair<double, double> interval(start, stop);
+	        std::vector < std::pair<double, double> > intervals;
+	        intervals.push_back(interval);
+	        frame_defs = stir::TimeFrameDefinitions(intervals);
+	        do_time_frame = true;
+        }
+
     private:
         //std::shared_ptr<PETAcquisitionData> sptr_ad_;
         std::shared_ptr<PETAcquisitionModelUsingMatrix> sptr_am_;

src/xSTIR/pSTIR/STIR.py

@@ -2333,6 +2333,15 @@ def set_up(self, image):
         """Sets up."""
         try_calling(pystir.cSTIR_setupPrior(self.handle, image.handle))
 
+    def accumulate_Hessian_times_input(self, current_estimate, input_, out=None):
+        """Computes the multiplication of the Hessian with a vector and adds it to output.
+        """
+        if out is None or out.handle is None:
+            out = input_.clone()
+        try_calling(pystir.cSTIR_priorAccumulateHessianTimesInput
+            (self.handle, out.handle, current_estimate.handle, input_.handle))
+        return out
+
 
 class QuadraticPrior(Prior):
     r"""Class for the prior that is a quadratic function of the image values.
@@ -2732,6 +2741,41 @@ def get_subset_gradient(self, image, subset, out=None):
         """
         return self.gradient(image, subset, out)
 
+    def accumulate_Hessian_times_input(self, current_estimate, input_, subset=-1, out=None):
+        """Computes the multiplication of the Hessian with a vector and adds it to output.
+        """
+        if out is None or out.handle is None:
+            out = input_.clone()
+        try_calling(pystir.cSTIR_objectiveFunctionAccumulateHessianTimesInput
+            (self.handle, current_estimate.handle, input_.handle, subset, out.handle))
+        return out
+
+    def multiply_with_Hessian(self, current_estimate, input_, subset=-1, out=None):
+        """Computes the multiplication of the Hessian with a vector.
+        """
+        if out is None or out.handle is None:
+            out = input_.clone()
+        try_calling(pystir.cSTIR_objectiveFunctionComputeHessianTimesInput
+            (self.handle, current_estimate.handle, input_.handle, subset, out.handle))
+        return out
+
+    def test_Hessian(self, x, subset=-1, eps=1e-3):
+        """Checks that grad(x + dx) - grad(x) is close to H(x)*dx
+        """
+        dx = x.clone()
+        dx *= eps/dx.norm()
+        dx += eps/2
+        y = x + dx
+        gx = self.gradient(x, subset)
+        gy = self.gradient(y, subset)
+        dg = gy - gx
+        Hdx = self.multiply_with_Hessian(x, dx, subset)
+        q = 1 - dg.norm()/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: %f' % q)
+        return q
+
     @abc.abstractmethod
     def get_subset_sensitivity(self, subset):
         #print('in base class ObjectiveFunction')
@@ -2854,7 +2898,8 @@ def set_acquisition_data(self, ad):
             self.handle, self.name, 'acquisition_data', ad.handle)
 
 
-class PoissonLogLikelihoodWithLinearModelForMeanAndListModeDataWithProjMatrixByBin(ObjectiveFunction):
+class PoissonLogLikelihoodWithLinearModelForMeanAndListModeDataWithProjMatrixByBin(PoissonLogLikelihoodWithLinearModelForMean):
+#(ObjectiveFunction):
     """Class for a STIR type of Poisson loglikelihood object for listmode data.
 
     Specifically, PoissonLogLikelihoodWithLinearModelForMeanAndListModeDataWithProjMatrixByBin.
@@ -2878,6 +2923,17 @@ def set_cache_path(self, path):
     def get_cache_path(self):
         return parms.char_par(self.handle, self.name, 'cache_path')
 
+    def set_time_interval(self, start, stop):
+        """Sets the time interval.
+
+        Only data scanned during this time interval will be converted.
+        """
+        interval = numpy.ndarray((2,), dtype=numpy.float32)
+        interval[0] = start
+        interval[1] = stop
+        try_calling(pystir.cSTIR_objFunListModeSetInterval(
+            self.handle, interval.ctypes.data))
+
     def set_acquisition_data(self, ad):
         assert_validity(ad, ListmodeData)
         parms.set_parameter(
@@ -2918,6 +2974,9 @@ def set_cache_max_size(self, diff):
     def get_cache_max_size(self):
         return parms.int_par(self.handle, self.name, 'cache_max_size')
 
+    def set_subsensitivity_filenames(self, names):
+        return parms.set_char_par(self.handle, self.name, 'subsensitivity_filenames', names)
+
     def get_subsensitivity_filenames(self):
         return parms.char_par(self.handle, self.name, 'subsensitivity_filenames')