Development branch for SQN research (DO NOT MERGE)

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/GeneralisedPrior.h

@@ -60,6 +60,17 @@ class GeneralisedPrior:
   virtual void compute_gradient(DataT& prior_gradient, 
 		   const DataT &current_estimate) =0; 
 
+  //! This should computes a single row of the Hessian
+  /*! Default implementation just call error(). This function needs to be overridden by the
+      derived class.
+      This method computes the Hessian of a particular voxel, indicated by \c coords, of the current image estimate.
+      The default method should overwrite the values in \c prior_Hessian_for_single_densel.
+   */
+  virtual Succeeded
+    compute_Hessian(DataT& prior_Hessian_for_single_densel,
+                    const BasicCoordinate<3,int>& coords,
+                    const DataT& current_image_estimate) const;
+
   //! This should compute the multiplication of the Hessian with a vector and add it to \a output
   /*! Default implementation just call error(). This function needs to be overridden by the
       derived class.
@@ -89,6 +100,9 @@ class GeneralisedPrior:
   virtual Succeeded 
     set_up(shared_ptr<const DataT> const& target_sptr);
 
+  //! Returns the status of the _is_convex variable
+  bool get_is_convex() const;
+
 protected:
   float penalisation_factor;
   //! sets value for penalisation factor
@@ -102,6 +116,9 @@ class GeneralisedPrior:
   virtual void check(DataT const& current_estimate) const;
 
   bool _already_set_up;
+
+  //! Variable to indicate that the prior is a convex function, should be set in defaults by the derived class
+  bool _is_convex;
 };
 
 END_NAMESPACE_STIR

src/include/stir/recon_buildblock/LogcoshPrior.h

@@ -125,9 +125,10 @@ class LogcoshPrior:  public
                                        const DiscretisedDensity<3,elemT> &current_image_estimate);
 
     //! compute Hessian
-    void compute_Hessian(DiscretisedDensity<3,elemT>& prior_Hessian_for_single_densel,
-                         const BasicCoordinate<3,int>& coords,
-                         const DiscretisedDensity<3,elemT> &current_image_estimate);
+    virtual Succeeded
+    compute_Hessian(DiscretisedDensity<3,elemT>& prior_Hessian_for_single_densel,
+                    const BasicCoordinate<3,int>& coords,
+                    const DiscretisedDensity<3,elemT> &current_image_estimate) const;
 
     //! Compute the multiplication of the hessian of the prior multiplied by the input.
     virtual Succeeded accumulate_Hessian_times_input(DiscretisedDensity<3,elemT>& output,
@@ -224,18 +225,19 @@ class LogcoshPrior:  public
       { return tanh(x)/x; }
     }
 
-    //! The Hessian of log(cosh()) is sech^2(x) = (1/cosh(x))^2
+    //! The second partial derivatives of the LogCosh Prior
     /*!
-     This function returns the hessian of the logcosh function
-     * @param d the difference between the ith and jth voxel.
-     * @param scalar is the logcosh scalar value controlling the priors transition between the quadratic and linear behaviour
-     * @return the second derivative of the log-cosh function
+     Diagonal refers to the second derivative w.r.t. x_j only (i.e. diagonal of the Hessian)
+     Off-diagonal refers to the second derivative w.r.t. x_j and x_k (i.e. off-diagonal of the Hessian)
+     For LogCosh, the off diagonal is the negative of the diagonal.
+     * @param x_j is the target voxel.
+     * @param x_k is the voxel in the neighbourhood.
+     * @return the second order partial derivatives of the LogCosh Prior
      */
-    static inline float Hessian(const float d, const float scalar)
-    {
-      const float x = d * scalar;
-      return square((1/ cosh(x)));
-    }
+    //@{
+    float diagonal_second_derivative(const float x_j, const float x_k) const;
+    float off_diagonal_second_derivative(const float x_j, const float x_k) const;
+    //@}
 };
 
 

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,