Merge pull request #1340 from KrisThielemans/already_set_up_obj_fun

Add already_set_up to objective function hierarchy

documentation/release_6.0.htm

@@ -252,6 +252,11 @@ <h3>Backward incompatibities</h3>
         <code>ProjDataInfo*NoArcCorr::get_bin_for_det_pair</code> is now private.
         Use <code>get_bin_for_det_pos_pair</code> instead.
       </li>
+      <li>
+        The <code>GeneralisedObjectiveFunction</code> hierarchy now has a <code>already_set_up</code>
+        member variable that needs to be set to <code>false</code> by <code>set_*</code>
+        functions and checked by callers.
+      </li>
     </ul>
 
     <h3>New functionality</h3>

recon_test_pack/run_tests.sh

@@ -4,7 +4,7 @@
 #  Copyright (C) 2000 - 2001 PARAPET partners
 #  Copyright (C) 2001 - 2009-10-11, Hammersmith Imanet Ltd
 #  Copyright (C) 2011, Kris Thielemans
-#  Copyright (C) 2013 - 2014, University College London
+#  Copyright (C) 2013 - 2014, 2024, University College London
 #  This file is part of STIR.
 #
 #  SPDX-License-Identifier: Apache-2.0 AND License-ref-PARAPET-license
@@ -19,7 +19,7 @@ if [ -n "$TRAVIS" -o -n "$GITHUB_WORKSPACE" ]; then
     set -e
 fi
 
-echo This script should work with STIR version 5.2. If you have
+echo This script should work with STIR version 6.0. If you have
 echo a later version, you might have to update your test pack.
 echo Please check the web site.
 echo
@@ -96,16 +96,20 @@ echo --------- TESTS THAT USE INTERPOLATING BACKPROJECTOR --------
 echo
 echo ------------- Running OSMAPOSL for sensitivity ------------- 
 echo Running ${INSTALL_DIR}OSMAPOSL for sensitivity
-${MPIRUN} ${INSTALL_DIR}OSMAPOSL OSMAPOSL_test_for_sensitivity.par 1> OSMAPOSL_test_for_sensitivity.log 2> OSMAPOSL_test_for_sensitivity_stderr.log 
-
-echo '---- Comparing output of sensitivity (should be identical up to tolerance)'
-echo Running ${INSTALL_DIR}compare_image
-if ${INSTALL_DIR}compare_image RPTsens_seg4.hv my_RPTsens_seg4.hv;
+if ${MPIRUN} ${INSTALL_DIR}OSMAPOSL OSMAPOSL_test_for_sensitivity.par 1> OSMAPOSL_test_for_sensitivity.log 2> OSMAPOSL_test_for_sensitivity_stderr.log
 then
-echo ---- This test seems to be ok !;
+    echo '---- Comparing output of sensitivity (should be identical up to tolerance)'
+    echo Running ${INSTALL_DIR}compare_image
+    if ${INSTALL_DIR}compare_image RPTsens_seg4.hv my_RPTsens_seg4.hv;
+    then
+        echo ---- This test seems to be ok !;
+    else
+        echo There were problems here!;
+        ThereWereErrors=1;
+    fi
 else
-echo There were problems here!;
-ThereWereErrors=1;
+    echo There were problems here!;
+    ThereWereErrors=1;
 fi
 
 echo
@@ -143,18 +147,23 @@ ${INSTALL_DIR}generate_image generate_uniform_image.par
 ${INSTALL_DIR}postfilter my_uniform_image_circular.hv my_uniform_image.hv postfilter_truncate_circular_FOV.par
 echo ------------- Running OSMAPOSL for sensitivity ------------- 
 echo Running ${INSTALL_DIR}OSMAPOSL for sensitivity
-${MPIRUN} ${INSTALL_DIR}OSMAPOSL OSMAPOSL_test_PM_for_sensitivity.par 1> sensitivity_PM.log 2> sensitivity_PM_stderr.log
-
-echo '---- Comparing output of sensitivity (should be identical up to tolerance)'
-echo Running ${INSTALL_DIR}compare_image
-if ${INSTALL_DIR}compare_image RPTsens_seg3_PM.hv my_RPTsens_seg3_PM.hv;
+if ${MPIRUN} ${INSTALL_DIR}OSMAPOSL OSMAPOSL_test_PM_for_sensitivity.par 1> sensitivity_PM.log 2> sensitivity_PM_stderr.log
 then
-echo ---- This test seems to be ok !;
+    echo '---- Comparing output of sensitivity (should be identical up to tolerance)'
+    echo Running ${INSTALL_DIR}compare_image
+    if ${INSTALL_DIR}compare_image RPTsens_seg3_PM.hv my_RPTsens_seg3_PM.hv;
+    then
+        echo ---- This test seems to be ok !;
+    else
+        echo There were problems here!;
+        ThereWereErrors=1;
+    fi
 else
-echo There were problems here!;
-ThereWereErrors=1;
+    echo There were problems here!;
+    ThereWereErrors=1;
 fi
 
+
 echo
 echo -------- Running OSMAPOSL  with the MRP prior -------- 
 echo Running ${INSTALL_DIR}OSMAPOSL

src/buildblock/ProjDataInfoGenericNoArcCorr.cxx

@@ -74,18 +74,8 @@ ProjDataInfoGenericNoArcCorr(const shared_ptr<Scanner> scanner_sptr,
   this->initialise_det1det2_to_uncompressed_view_tangpos();
 #endif
 
-  CartesianCoordinate3D< float> b1,b2;
-  Bin bin;
-  bin.segment_num() = 0;
-  bin.axial_pos_num() = 0;
-  bin.view_num() = 0;
-  bin.tangential_pos_num() = 0;
-// setting shift_z to 0 before it is actually estimated. Otherwise the next function will use it
-  this->z_shift.z()=0;
-  find_cartesian_coordinates_of_detection(b1,b2,bin);
-  float shift=b2.z();
-
-  this->z_shift.z()=shift;
+  // find shift between "new" centre-of-scanner and "old" centre-of-first-ring coordinate system
+  this->z_shift.z()=this->get_scanner_ptr()->get_coordinate_for_det_pos(DetectionPosition<>(0,0,0)).z();
   this->z_shift.y()=0;
   this->z_shift.x()=0;
 }

src/experimental/motion/PoissonLogLikelihoodWithLinearModelForMeanAndGatedProjDataWithMotion.cxx

@@ -285,6 +285,7 @@ void
 PoissonLogLikelihoodWithLinearModelForMeanAndGatedProjDataWithMotion<TargetT>::
 set_proj_data_sptr(const shared_ptr<GatedProjData>& arg)
 {
+  this->already_set_up = false;
   this->_gated_proj_data_sptr = arg;
 }
 
@@ -293,15 +294,16 @@ void
 PoissonLogLikelihoodWithLinearModelForMeanAndGatedProjDataWithMotion<TargetT>::
 set_max_segment_num_to_process(const int arg)
 {
+  this->already_set_up = this->already_set_up && (this->max_timing_pos_num_to_process == arg);
   this->max_segment_num_to_process = arg;
-
 }
 
 template<typename TargetT>
 void
 PoissonLogLikelihoodWithLinearModelForMeanAndGatedProjDataWithMotion<TargetT>::
 set_zero_seg0_end_planes(const bool arg)
 {
+  this->already_set_up = this->already_set_up && (this->zero_seg0_end_planes == arg);
   this->zero_seg0_end_planes = arg;
 }
 
@@ -310,7 +312,7 @@ void
 PoissonLogLikelihoodWithLinearModelForMeanAndGatedProjDataWithMotion<TargetT>::
 set_additive_proj_data_sptr(const shared_ptr<GatedProjData>& arg)
 {
-
+  this->already_set_up = false;
   this->_gated_additive_proj_data_sptr = arg;
 }
 
@@ -319,6 +321,7 @@ void
 PoissonLogLikelihoodWithLinearModelForMeanAndGatedProjDataWithMotion<TargetT>::
 set_projector_pair_sptr(const shared_ptr<ProjectorByBinPair>& arg) 
 {
+  this->already_set_up = false;  
   this->projector_pair_ptr = arg;
 
 }
@@ -328,6 +331,7 @@ void
 PoissonLogLikelihoodWithLinearModelForMeanAndGatedProjDataWithMotion<TargetT>::
 set_frame_num(const int arg)
 {
+  this->already_set_up = this->already_set_up && (this->frame_num == arg);
   this->frame_num = arg;
 }
 
@@ -336,6 +340,7 @@ void
 PoissonLogLikelihoodWithLinearModelForMeanAndGatedProjDataWithMotion<TargetT>::
 set_frame_definitions(const TimeFrameDefinitions& arg)
 {
+  this->already_set_up = this->already_set_up && (this->frame_defs == arg);
   this->frame_defs = arg;
 }
 

src/include/stir/recon_buildblock/GeneralisedObjectiveFunction.h

@@ -86,12 +86,19 @@ class GeneralisedObjectiveFunction:
 {
 public:
 
-  //GeneralisedObjectiveFunction(); 
+  GeneralisedObjectiveFunction()
+    : already_set_up(false)
+  {}
+
 
   virtual ~GeneralisedObjectiveFunction(); 
 
 
   //! Creates a suitable target as determined by the parameters
+  /*!
+    \warning This should <b>not</b> check \c already_set_up (unfortunately),
+    as it is currently called in Reconstruction::reconstruct() before calling set_up().
+  */
   virtual TargetT *
     construct_target_ptr() const = 0; 
 
@@ -331,6 +338,7 @@ class GeneralisedObjectiveFunction:
 
 protected:
   int num_subsets;
+  bool already_set_up;
 
   shared_ptr<GeneralisedPrior<TargetT> > prior_sptr;
 

src/include/stir/recon_buildblock/PoissonLogLikelihoodWithLinearModelForMeanAndGatedProjDataWithMotion.txx

@@ -175,7 +175,7 @@ template <typename TargetT>
 TargetT *
 PoissonLogLikelihoodWithLinearModelForMeanAndGatedProjDataWithMotion<TargetT>::
 construct_target_ptr() const
-{  
+{
   return 
     this->target_parameter_parser.create(this->get_input_data());
 }
@@ -280,29 +280,34 @@ get_projector_pair_sptr() const
 template<typename TargetT>
 int
 PoissonLogLikelihoodWithLinearModelForMeanAndGatedProjDataWithMotion<TargetT>::
-set_num_subsets(const int num_subsets)
+set_num_subsets(const int new_num_subsets)
 {
+  this->already_set_up = this->already_set_up && (this->num_subsets == new_num_subsets);
   for(unsigned int gate_num=1;gate_num<=this->get_time_gate_definitions().get_num_gates();++gate_num)
     {
       if(this->_single_gate_obj_funcs.size() != 0)
-		if(this->_single_gate_obj_funcs[gate_num].set_num_subsets(num_subsets) != num_subsets)
+		if(this->_single_gate_obj_funcs[gate_num].set_num_subsets(new_num_subsets) != new_num_subsets)
 			error("set_num_subsets didn't work");
     }
-  this->num_subsets=num_subsets;
+  this->num_subsets=new_num_subsets;
   return this->num_subsets;
 }
 
 template<typename TargetT>
 void
 PoissonLogLikelihoodWithLinearModelForMeanAndGatedProjDataWithMotion<TargetT>::
 set_time_gate_definitions(const TimeGateDefinitions & time_gate_definitions)
-{ this->_time_gate_definitions=time_gate_definitions; }
+{
+  this->already_set_up = false;
+  this->_time_gate_definitions=time_gate_definitions;
+}
 
 template<typename TargetT>
 void
 PoissonLogLikelihoodWithLinearModelForMeanAndGatedProjDataWithMotion<TargetT>::
 set_input_data(const shared_ptr<ExamData> & arg)
 {
+  this->already_set_up = false;
     this->_gated_proj_data_sptr = dynamic_pointer_cast<GatedProjData>(arg);
 }
 
@@ -319,6 +324,7 @@ void
 PoissonLogLikelihoodWithLinearModelForMeanAndGatedProjDataWithMotion<TargetT>::
 set_additive_proj_data_sptr(const shared_ptr<ExamData> &arg)
 {
+  this->already_set_up = false;
     this->_additive_gated_proj_data_sptr = dynamic_pointer_cast<GatedProjData>(arg);
 }
 
@@ -327,6 +333,7 @@ void
 PoissonLogLikelihoodWithLinearModelForMeanAndGatedProjDataWithMotion<TargetT>::
 set_normalisation_sptr(const shared_ptr<BinNormalisation>& arg)
 {
+  this->already_set_up = false;
 //  this->normalisation_sptr = arg;
     error("Not implemeted yet");
 }

src/recon_buildblock/GeneralisedObjectiveFunction.cxx

@@ -142,6 +142,8 @@ compute_sub_gradient(TargetT& gradient,
 		     const TargetT &current_estimate, 
 		     const int subset_num)
 {
+  if (!this->already_set_up)
+    error("Need to call set_up() for objective function first");
   assert(gradient.get_index_range() == current_estimate.get_index_range());
 
   if (subset_num<0 || subset_num>=this->get_num_subsets())
@@ -194,6 +196,8 @@ GeneralisedObjectiveFunction<TargetT>::
 compute_objective_function_without_penalty(const TargetT& current_estimate,
 					   const int subset_num)
 {
+  if (!this->already_set_up)
+    error("Need to call set_up() for objective function first");
   if (subset_num<0 || subset_num>=this->get_num_subsets())
     error("compute_objective_function_without_penalty subset_num out-of-range error");
 
@@ -231,6 +235,8 @@ add_multiplication_with_approximate_sub_Hessian_without_penalty(TargetT& output,
 								const TargetT& input,
 								const int subset_num) const
 {
+  if (!this->already_set_up)
+    error("Need to call set_up() for objective function first");
   if (subset_num<0 || subset_num>=this->get_num_subsets())
     error("add_multiplication_with_approximate_sub_Hessian_without_penalty subset_num out-of-range error");
 
@@ -259,6 +265,8 @@ add_multiplication_with_approximate_sub_Hessian(TargetT& output,
 						const TargetT& input,
 						const int subset_num) const
 {
+  if (!this->already_set_up)
+    error("Need to call set_up() for objective function first");
   if (this->add_multiplication_with_approximate_sub_Hessian_without_penalty(output, input, subset_num) ==
       Succeeded::no)
     return Succeeded::no;
@@ -406,6 +414,8 @@ accumulate_sub_Hessian_times_input_without_penalty(TargetT& output,
                                                    const TargetT& input,
                                                    const int subset_num) const
 {
+  if (!this->already_set_up)
+    error("Need to call set_up() for objective function first");
   if (subset_num<0 || subset_num>=this->get_num_subsets())
     error("accumulate_sub_Hessian_times_input_without_penalty subset_num out-of-range error");
 
@@ -487,6 +497,12 @@ bool
 GeneralisedObjectiveFunction<TargetT>::
 subsets_are_approximately_balanced(std::string& warning_message) const
 {
+#if 0
+  // TODO cannot do this yet, as this function is called in
+  // `PoissonLogLikelihoodWithLinearModelForMean::compute_sensitivities` during set_up()
+  if (!this->already_set_up)
+    error("Need to call set_up() for objective function first");
+#endif
   return this->actual_subsets_are_approximately_balanced(warning_message);
 }
 

src/recon_buildblock/PoissonLogLikelihoodWithLinearModelForMean.cxx

@@ -95,6 +95,7 @@ void
 PoissonLogLikelihoodWithLinearModelForMean<TargetT>::
 set_sensitivity_filename(const std::string& filename)
 {
+  this->already_set_up = false;
   this->sensitivity_filename = filename;
 }
 
@@ -103,6 +104,7 @@ void
 PoissonLogLikelihoodWithLinearModelForMean<TargetT>::
 set_subsensitivity_filenames(const std::string& filenames)
 {
+  this->already_set_up = false;
   this->subsensitivity_filenames = filenames;
   try
     {
@@ -171,6 +173,7 @@ void
 PoissonLogLikelihoodWithLinearModelForMean<TargetT>::
 set_use_subset_sensitivities(const bool arg)
 {
+  this->already_set_up = this->already_set_up && (this->use_subset_sensitivities == arg);
   this->use_subset_sensitivities = arg;
 }
 
@@ -179,6 +182,7 @@ void
 PoissonLogLikelihoodWithLinearModelForMean<TargetT>::
 set_subset_sensitivity_sptr(const shared_ptr<TargetT>& arg, const int subset_num)
 {
+  this->already_set_up = false;
   this->subsensitivity_sptrs[subset_num] = arg;
 }
 
@@ -348,7 +352,7 @@ set_up(shared_ptr<TargetT> const& target_sptr)
           return Succeeded::no;
         }
     }
-      
+  this->already_set_up = true;
   return Succeeded::yes;
 }
 
@@ -359,6 +363,8 @@ compute_sub_gradient_without_penalty(TargetT& gradient,
                                      const TargetT &current_estimate, 
                                      const int subset_num)
 {
+  if (!this->already_set_up)
+    error("Need to call set_up() for objective function first");
   this->actual_compute_subset_gradient_without_penalty(gradient, current_estimate, subset_num, false);
 }
 
@@ -369,6 +375,8 @@ compute_sub_gradient_without_penalty_plus_sensitivity(TargetT& gradient,
                                                       const TargetT &current_estimate,
                                                       const int subset_num)
 {
+  if (!this->already_set_up)
+    error("Need to call set_up() for objective function first");
   actual_compute_subset_gradient_without_penalty(gradient, current_estimate, subset_num, true);
 }