Removing requirement for gradient to use sensitivity and instead does…

… proper gradient computation (#19)

* Create a templated method of RPC_process_related_viewgrams_gradient

* Update the standard Likelihood classes to use actual_compute_sub_gradient_without_penalty

* Dynamic data actual_compute_sub_gradient_without_penalty

* Gated data actual_compute_sub_gradient_without_penalty

* List mode actual_compute_sub_gradient_without_penalty will error if do_subtraction is true

* Remove overrides [ci skip]

* Update documentation

* Gradient computation will subtract normalised ones (mult_viewgrams)

* [ci skip] Minor commit

* If zero_seg0_end_planes and no normalisation, create mult_viewgrams so the endplanes can be zeroed

* remove compute_sub_gradient_without_penalty_plus_sensitivity from child Poisson classes

* Fix documentation formatting and code whitespace

* Gated and Dynamic actual's should call the objective function actual

* reformat do_subtraction to add_sensitivity (involves inverting cases)

* Revert one too many inverts

* Improve comments

* List mode recon: !add_sensitivity and error if not using subset sensitivity

* turn off add_sensitivity on compute_sub_gradient_without_penalty

* distributable_compute_gradient: swap normalisation_sptr parsing case

* Refactor actual_compute_sub_grad... to actual_compute_subset_grad...

* Update OSSPS test image correct OSSPS (with actual subset sensitivity subtraction)

* Add OSSPS subset sensitivity test. Results should be identical

* [ci skip] Add compute_sub_gradient_without_penalty documentation

* [ci skip] Documentation cleanup

* [ci skip] Add to release notes

documentation/release_5.0.htm

@@ -143,6 +143,12 @@ <h3>Changed functionality</h3>
 	<li><tt>ROOT</tt> file I/O improvements. An entire entry's tree is no longer loaded by default and instead individual branches are loaded as needed.
     ROOT file event processing is now up to 4x faster. In addition, there is now an option to <tt>check energy window information</tt> (defaulting to on).
     Futhermore, added functionality to exclude true and unscattered event types from list mode processing. </li>
+    <li>
+        Poisson log-likelihood gradient methods now use <tt>actual_compute_subset_gradient_without_penalty</tt>
+        to properly handle the objective function gradient computation at the distributable level.
+        Removes the subtraction of the subset sensitivity, which lead to a inconsistency between log-likelihood function
+        and gradient when <tt>use_subset_sensitivities</tt> was false.
+    </li>
 </ul>
 
 

recon_test_pack/OSSPS_test_PM_QP_subsens1.par

@@ -0,0 +1,73 @@
+OSSPSParameters :=
+; sample file for OSSPS
+; parameters used here are for illustrative purposes only
+; i.e. they are not recommended values
+
+objective function type:= PoissonLogLikelihoodWithLinearModelForMeanAndProjData
+PoissonLogLikelihoodWithLinearModelForMeanAndProjData Parameters:=
+
+input file := Utahscat600k_ca_seg4.hs
+zero end planes of segment 0:= 1
+; if disabled, defaults to maximum segment number in the file
+maximum absolute segment number to process := 3
+
+; change to STIR 2.x default for compatibility 
+use subset sensitivities:=1
+sensitivity filename:= RPTsens_seg3_PM.hv
+
+projector pair type := Matrix
+  Projector Pair Using Matrix Parameters :=
+  Matrix type := Ray Tracing
+  Ray tracing matrix parameters :=
+   number of rays in tangential direction to trace for each bin := 2
+;   restrict to cylindrical fov := 0
+  End Ray tracing matrix parameters :=
+  End Projector Pair Using Matrix Parameters :=
+
+; additive sinogram:=my_fake_randoms.hs
+prior type := quadratic
+  Quadratic Prior Parameters:=
+  penalisation factor := 0.5
+  ; next defaults to 0, set to 1 for 2D inverse Euclidean weights, 0 for 3D 
+  only 2D:= 0
+  END Quadratic Prior Parameters:=
+
+End PoissonLogLikelihoodWithLinearModelForMeanAndProjData Parameters:=
+
+output filename prefix := my_test_image_OSSPS_PM_QP_subsens
+; iteration scheme
+
+number of subsets:= 4
+;start at subset:= 0
+;start at subiteration number := 1
+number of subiterations:= 8
+Save estimates at subiteration intervals:= 8
+;write update image := 0
+;report objective function values interval := 1
+
+; if next is disabled, defaults to image full of 1s (but that's not good for OSSPS)
+; in particular, make sure it has the correct scale
+initial estimate:= test_image_PM_QP_6.hv
+enforce initial positivity condition := 1
+
+; here start OSSPS specific values
+
+; values to use for the 'precomputed denominator'
+; specify either procomputed denomiator or normalisation type
+  ; use the following if you have it already (e.g. from previous run)
+  ; note: setting the value to 1 will use an images full of ones.
+  ; precomputed denominator := my_precomputed_denominator.hv
+
+; specify relaxation scheme
+; lambda = relaxation_parameter/ (1+relaxation_gamma*(subiteration_num/num_subsets)
+relaxation parameter := 2
+relaxation gamma:=.1
+
+
+END :=
+
+
+
+
+
+

recon_test_pack/run_tests.sh

@@ -245,6 +245,21 @@ echo There were problems here!;
 ThereWereErrors=1;
 fi
 
+${MPIRUN} ${INSTALL_DIR}OSSPS OSSPS_test_PM_QP_subsens1.par 1> OSSPS_PM_QP.log 2> OSSPS_PM_QP_stderr.log
+
+echo '---- Comparing output of OSSPS subiter 8 using subset sensitivity (should be identical up to tolerance)'
+echo Running ${INSTALL_DIR}compare_image
+# relax test for the outer-rim voxels as these turn out to be more unstable than the internal ones
+if ${INSTALL_DIR}compare_image -t 0.002 test_image_OSSPS_PM_QP_8.hv my_test_image_OSSPS_PM_QP_subsens_8.hv -a
+   ${INSTALL_DIR}compare_image -r 1 test_image_OSSPS_PM_QP_8.hv my_test_image_OSSPS_PM_QP_subsens_8.hv
+then
+  echo ---- This test seems to be ok !;
+  else
+  echo There were problems here!;
+  ThereWereErrors=1;
+fi
+
+
 echo
 echo ------------- tests on stir_math and correct_projdata ---------
   echo "first make up some randoms (just a projdata full of 1)"

src/include/stir/recon_buildblock/PoissonLogLikelihoodWithLinearKineticModelAndDynamicProjectionData.h

@@ -70,10 +70,6 @@ public  RegisteredParsingObject<PoissonLogLikelihoodWithLinearKineticModelAndDyn
   */
   virtual TargetT *
     construct_target_ptr() const; 
-  virtual void 
-    compute_sub_gradient_without_penalty_plus_sensitivity(TargetT& gradient, 
-                                                          const TargetT &current_estimate, 
-                                                          const int subset_num); 
 
   virtual std::unique_ptr<ExamInfo>
   get_exam_info_uptr_for_target()  const;
@@ -170,6 +166,13 @@ public  RegisteredParsingObject<PoissonLogLikelihoodWithLinearKineticModelAndDyn
 
   bool actual_subsets_are_approximately_balanced(std::string& warning_message) const;
 
+
+  virtual void
+  actual_compute_subset_gradient_without_penalty(TargetT& gradient,
+                                                 const TargetT &current_estimate,
+                                                 const int subset_num,
+                                                 const bool add_sensitivity);
+
   //! Sets defaults for parsing 
   /*! Resets \c sensitivity_filename and \c sensitivity_sptr and
      \c recompute_sensitivity to \c false.

src/include/stir/recon_buildblock/PoissonLogLikelihoodWithLinearKineticModelAndDynamicProjectionData.txx

@@ -426,9 +426,10 @@ set_normalisation_sptr(const shared_ptr<BinNormalisation>& arg)
 template<typename TargetT>
 void
 PoissonLogLikelihoodWithLinearKineticModelAndDynamicProjectionData<TargetT>::
-compute_sub_gradient_without_penalty_plus_sensitivity(TargetT& gradient, 
-                                                      const TargetT &current_estimate, 
-                                                      const int subset_num)
+actual_compute_subset_gradient_without_penalty(TargetT& gradient,
+                                               const TargetT &current_estimate,
+                                               const int subset_num,
+                                               const bool add_sensitivity)
 {
   if (subset_num<0 || subset_num>=this->get_num_subsets())
     error("compute_sub_gradient_without_penalty subset_num out-of-range error");
@@ -452,9 +453,10 @@ compute_sub_gradient_without_penalty_plus_sensitivity(TargetT& gradient,
 
 
       this->_single_frame_obj_funcs[frame_num].
-        compute_sub_gradient_without_penalty_plus_sensitivity(dyn_gradient[frame_num], 
-                                                              dyn_image_estimate[frame_num], 
-                                                              subset_num);
+              actual_compute_subset_gradient_without_penalty(dyn_gradient[frame_num],
+                                                             dyn_image_estimate[frame_num],
+                                                             subset_num,
+                                                             add_sensitivity);
     }
 
   this->_patlak_plot_sptr->multiply_dynamic_image_with_model_gradient(gradient,

src/include/stir/recon_buildblock/PoissonLogLikelihoodWithLinearModelForMean.h

@@ -66,14 +66,13 @@ START_NAMESPACE_STIR
   the <i>sensitivity</i> because (if \f$r=0\f$) it is the total
   probability of detecting a count (in any bin) originating from \f$v\f$.
 
-  This class computes the gradient as a sum of these two terms. The
-  sensitivity has to be computed by the virtual function 
-  \c add_subset_sensitivity(). The sum is computed by
-  \c compute_sub_gradient_without_penalty_plus_sensitivity().
-
-  The reason for this is that the sensitivity is data-independent, and
-  can be computed only once. See also
-  PoissonLogLikelihoodWithLinearModelForMeanAndListModeData.
+  This class computes the gradient directly, via \c compute_sub_gradient_without_penalty().
+  This method is utilised by the \c OSSPS algorithm in STIR.
+  However, an additional method (\c compute_sub_gradient_without_penalty_plus_sensitivity())
+  is provided that computes the sum of the subset gradient (without penalty) and the sensitivity.
+  This method is utilised by the \c OSMAPOSL algorithm.
+
+  See also \c PoissonLogLikelihoodWithLinearModelForMeanAndListModeData.
   
   \par Relation with Kullback-Leibler distance
 
@@ -120,38 +119,40 @@ public  GeneralisedObjectiveFunction<TargetT>
 
   //PoissonLogLikelihoodWithLinearModelForMean(); 
 
-  //! Implementation in terms of compute_sub_gradient_without_penalty_plus_sensitivity()
-  /*! \warning If separate subsensitivities are not used, we just subtract the total 
-    sensitivity divided by the number of subsets.
-    This is fine for some algorithms as the sum over all the subsets is 
-    equal to gradient of the objective function (without prior). 
-    Other algorithms do not behave very stable under this approximation
-    however. So, currently setup() will return an error if
-    <code>!subsets_are_approximately_balanced()</code> and subset sensitivities
-    are not used.
-
-    \see get_use_subset_sensitivities()
-  */
+    //! Compute the subset gradient of the (unregularised) objective function
+    /*!
+     Implementation in terms of actual_compute_sub_gradient_without_penalty()
+     This function is used by OSSPS may be used by other gradient ascent/descent algorithms
+
+      This computes
+      \f[
+      {\partial L \over \partial \lambda_v} =
+        \sum_b P_{bv} ({y_b \over Y_b} - 1)
+      \f]
+      (see the class general documentation).
+      The sum will however be restricted to a subset.
+    */
   virtual void 
     compute_sub_gradient_without_penalty(TargetT& gradient, 
                                          const TargetT &current_estimate, 
                                          const int subset_num); 
 
-  //! This should compute the gradient of the (unregularised) objective function plus the (sub)sensitivity
-  /*! 
-    This function is used for instance by OSMAPOSL.
-
-    This computes
-    \f[ {\partial L \over \partial \lambda_v} + P_v =
-      \sum_b P_{bv} {y_b \over Y_b}
-      \f]
-    (see the class general documentation).
-    The sum will however be restricted to a subset.
-   */
+    //! This should compute the subset gradient of the (unregularised) objective function plus the subset sensitivity
+    /*!
+     Implementation in terms of actual_compute_sub_gradient_without_penalty().
+     This function is used for instance by OSMAPOSL.
+
+      This computes
+      \f[ {\partial L \over \partial \lambda_v} + P_v =
+        \sum_b P_{bv} {y_b \over Y_b}
+        \f]
+      (see the class general documentation).
+      The sum will however be restricted to a subset.
+     */
   virtual void 
     compute_sub_gradient_without_penalty_plus_sensitivity(TargetT& gradient, 
                                                           const TargetT &current_estimate, 
-                                                          const int subset_num) =0; 
+                                                          const int subset_num);
 
   //! set-up sensitivity etc if possible
   /*! If \c recompute_sensitivity is \c false, we will try to
@@ -264,6 +265,27 @@ public  GeneralisedObjectiveFunction<TargetT>
   */
   void compute_sensitivities();
 
+  //! computes the subset gradient of the objective function without the penalty (optional: add subset sensitivity)
+  /*!
+    If \c add_sensitivity is \c true, this computes
+    \f[ {\partial L \over \partial \lambda_v} + P_v =
+      \sum_b P_{bv} {y_b \over Y_b}
+      \f]
+    (see the class general documentation).
+    The sum will however be restricted to a subset.
+
+    However, if \c add_sensitivity is \c false, this function will instead compute only the gradient
+    \f[
+        {\partial L \over \partial \lambda_v} =
+            \sum_b P_{bv} ({y_b \over Y_b} - 1)
+    \f]
+  */
+  virtual void
+  actual_compute_subset_gradient_without_penalty(TargetT& gradient,
+                                                 const TargetT &current_estimate,
+                                                 const int subset_num,
+                                                 const bool add_sensitivity) = 0;
+
   //! Sets defaults for parsing 
   /*! Resets \c sensitivity_filename, \c subset_sensitivity_filenames to empty,
      \c recompute_sensitivity to \c false, and \c use_subset_sensitivities to false.

src/include/stir/recon_buildblock/PoissonLogLikelihoodWithLinearModelForMeanAndGatedProjDataWithMotion.h

@@ -79,10 +79,11 @@ public  RegisteredParsingObject<PoissonLogLikelihoodWithLinearModelForMeanAndGat
   virtual TargetT *
     construct_target_ptr() const; 
 
-  virtual void 
-    compute_sub_gradient_without_penalty_plus_sensitivity(TargetT& gradient, 
-							  const TargetT &current_estimate, 
-							  const int subset_num); 
+  virtual void
+    actual_compute_subset_gradient_without_penalty(TargetT& gradient,
+                                                 const TargetT &current_estimate,
+                                                 const int subset_num,
+                                                 const bool add_sensitivity);
 
   virtual double
     actual_compute_objective_function_without_penalty(const TargetT& current_estimate,

src/include/stir/recon_buildblock/PoissonLogLikelihoodWithLinearModelForMeanAndGatedProjDataWithMotion.txx

@@ -415,9 +415,10 @@ set_up_before_sensitivity(shared_ptr<const TargetT > const& target_sptr)
 template<typename TargetT>
 void
 PoissonLogLikelihoodWithLinearModelForMeanAndGatedProjDataWithMotion<TargetT>::
-compute_sub_gradient_without_penalty_plus_sensitivity(TargetT& gradient, 
-						      const TargetT &current_estimate, 
-						      const int subset_num)
+actual_compute_subset_gradient_without_penalty(TargetT& gradient,
+                                               const TargetT &current_estimate,
+                                               const int subset_num,
+                                               const bool add_sensitivity)
 {
   assert(subset_num>=0);
   assert(subset_num<this->num_subsets);
@@ -435,9 +436,10 @@ compute_sub_gradient_without_penalty_plus_sensitivity(TargetT& gradient,
               gated_gradient[gate_num].end_all(),
               0.F);
     this->_single_gate_obj_funcs[gate_num].
-      compute_sub_gradient_without_penalty_plus_sensitivity(gated_gradient[gate_num], 
-                                                            gated_image_estimate[gate_num], 
-                                                            subset_num);
+            actual_compute_subset_gradient_without_penalty(gated_gradient[gate_num],
+                                                           gated_image_estimate[gate_num],
+                                                           subset_num,
+                                                           add_sensitivity);
   }	
   //	if(this->_motion_correction_type==-1)
   this->_reverse_motion_vectors.warp_image(gradient,gated_gradient) ; 

src/include/stir/recon_buildblock/PoissonLogLikelihoodWithLinearModelForMeanAndListModeData.h

@@ -78,10 +78,6 @@ public PoissonLogLikelihoodWithLinearModelForMean<TargetT>
   virtual Succeeded
    set_up(shared_ptr <TargetT > const& target_sptr);
 
-  virtual  
-  void compute_sub_gradient_without_penalty_plus_sensitivity(TargetT& gradient,  
-	    				 const TargetT &current_estimate,  
-					 const int subset_num)=0;  
   //! time frame definitions
   /*! \todo This is currently used to be able to compute the gradient for 
     one time frame. However, it probably does not belong here.

src/include/stir/recon_buildblock/PoissonLogLikelihoodWithLinearModelForMeanAndListModeDataWithProjMatrixByBin.h

@@ -66,11 +66,17 @@ typedef RegisteredParsingObject<PoissonLogLikelihoodWithLinearModelForMeanAndLis
 
   PoissonLogLikelihoodWithLinearModelForMeanAndListModeDataWithProjMatrixByBin<TargetT>(); 
 
-  //! This should compute the gradient of the objective function at the \a current_estimate overwriting \a gradient
-  virtual  
-  void compute_sub_gradient_without_penalty_plus_sensitivity(TargetT& gradient,  
-                                         const TargetT &current_estimate,  
-                                         const int subset_num);  
+  //! Computes the gradient of the objective function at the \a current_estimate overwriting \a gradient.
+  /*!
+   \warning If <code>add_sensitivity = false</code> and <code>use_subset_sensitivities = false</code> will return an error
+   because the gradient will not be correct. Try <code>use_subset_sensitivities = true</code>.
+   */
+    virtual
+    void actual_compute_subset_gradient_without_penalty(TargetT& gradient,
+                                                        const TargetT &current_estimate,
+                                                        const int subset_num,
+                                                        const bool add_sensitivity);
+
   virtual TargetT * construct_target_ptr() const;  
 
   int set_num_subsets(const int new_num_subsets);

src/include/stir/recon_buildblock/PoissonLogLikelihoodWithLinearModelForMeanAndProjData.h

@@ -212,10 +212,11 @@ public  RegisteredParsingObject<PoissonLogLikelihoodWithLinearModelForMeanAndPro
   virtual const ProjData& get_input_data() const;
 //@}
 
-  virtual void 
-    compute_sub_gradient_without_penalty_plus_sensitivity(TargetT& gradient, 
-                                                          const TargetT &current_estimate, 
-                                                          const int subset_num); 
+  virtual void
+  actual_compute_subset_gradient_without_penalty(TargetT& gradient,
+                                                 const TargetT &current_estimate,
+                                                 const int subset_num,
+                                                 const bool add_sensitivity);
 
   virtual std::unique_ptr<ExamInfo>
   get_exam_info_uptr_for_target()  const;