Revert "Merge pull request UCL#901 from robbietuk/RDP/HessianVectorPr…

…oduct"

This reverts commit 4d087f3, reversing
changes made to f98246f.

documentation/release_5.1.htm

@@ -25,8 +25,7 @@ <h2>Overall summary</h2>
 improvements to the documentation.
 </p>
 
-<p>This release contains mainly code written by Nikos Eftimiou (UPenn and MGH), with support by Kris Thielemans (UCL).
-  Additional changes to penalty functions were made by Robert Twyman (UCL).
+<p>This release contains mainly code written by Nikos Eftimiou (UPenn and MGH), with support by Kris Thielemans (UCL)
 </p>
 
 <h2>Patch release info</h2>
@@ -44,8 +43,6 @@ <h3>New functionality</h3>
   </li>    
   <li>Support for PENNPET Explorer listmode data (if proprietary libraries are found) by Nikos Efthimiou, see <a href="https://github.com/UCL/STIR/pull/1028/">PR #1028</a>. 
   </li>    
-  <li> Improvements to penalty functions including Hessian methods, see <a href="https://github.com/UCL/STIR/pull/901/">PR #901</a>.
-  </li>
 </ul>
 
 
@@ -69,11 +66,6 @@ <h3>Documentation changes</h3>
 <h3>recon_test_pack changes</h3>
 
 <h3>Other changes to tests</h3>
-  <li> Significant changes made to <code>test_priors</code> that tests the Hessian's convexity, given
-    <code>x(Hx) > 0</code>, and a perturbation response, using gradients, was added to determine the Hessian
-    (for a single densel) is within tolerance.
-    Tests the Quadratic, Relative Difference (in two configurations) and Log-Cosh Penalties (Robert Twyman, UCL).
-  </li>
 
 </body>
 

src/include/stir/recon_buildblock/FilterRootPrior.h

@@ -87,8 +87,6 @@ class FilterRootPrior:  public
   FilterRootPrior(shared_ptr<DataProcessor<DataT> >const&, 
 		  const float penalization_factor);
 
-  bool is_convex() const;
-
  //! compute the value of the function
   /*! \warning Generally there is no function associated to this prior,
     so we just return 0 and write a warning the first time it's called.

src/include/stir/recon_buildblock/GeneralisedPrior.h

@@ -60,27 +60,14 @@ class GeneralisedPrior:
   virtual void compute_gradient(DataT& prior_gradient, 
 		   const DataT &current_estimate) =0; 
 
-  //! This computes a single row of the Hessian
-  /*! Default implementation just call error(). This function needs to be overridden by the
-      derived class.
-      
-      The method computes a row (i.e. at a densel/voxel, indicated by \c coords) of the Hessian at \c current_estimate. 
-      Note that a row corresponds to an object of `DataT`.
-      The method (as implemented in derived classes) should store the result in \c prior_Hessian_for_single_densel.
-   */
-  virtual void
-    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.
       This method assumes that the hessian of the prior is 1 and hence the function quadratic.
       Instead, accumulate_Hessian_times_input() should be used. This method remains for backwards comparability.
        \warning The derived class should accumulate in \a output.
   */
-  virtual void
+  virtual Succeeded 
     add_multiplication_with_approximate_Hessian(DataT& output,
 						const DataT& input) const;
 
@@ -89,7 +76,7 @@ class GeneralisedPrior:
         derived class.
         \warning The derived class should accumulate in \a output.
     */
-  virtual void
+  virtual Succeeded
   accumulate_Hessian_times_input(DataT& output,
           const DataT& current_estimate,
           const DataT& input) const;
@@ -102,10 +89,6 @@ class GeneralisedPrior:
   virtual Succeeded 
     set_up(shared_ptr<const DataT> const& target_sptr);
 
-  //! Indicates if the prior is a smooth convex function
-  /*! If true, the prior is expected to have 0th, 1st and 2nd order behaviour implemented.*/
-  virtual bool is_convex() const = 0;
-
 protected:
   float penalisation_factor;
   //! sets value for penalisation factor

src/include/stir/recon_buildblock/LogcoshPrior.h

@@ -108,12 +108,6 @@ class LogcoshPrior:  public
     //! Constructs it explicitly with scalar
     LogcoshPrior(const bool only_2D, float penalization_factor, const float scalar);
 
-    virtual bool
-    parabolic_surrogate_curvature_depends_on_argument() const
-    { return false; }
-
-    bool is_convex() const;
-
     //! compute the value of the function
     double
     compute_value(const DiscretisedDensity<3,elemT> &current_image_estimate);
@@ -126,13 +120,13 @@ class LogcoshPrior:  public
     void parabolic_surrogate_curvature(DiscretisedDensity<3,elemT>& parabolic_surrogate_curvature,
                                        const DiscretisedDensity<3,elemT> &current_image_estimate);
 
-    virtual void
-    compute_Hessian(DiscretisedDensity<3,elemT>& prior_Hessian_for_single_densel,
-                    const BasicCoordinate<3,int>& coords,
-                    const DiscretisedDensity<3,elemT> &current_image_estimate) const;
+    //! 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);
 
-    //! Compute the multiplication of the hessian of the prior (at \c current_estimate) and the \c input.
-    virtual void accumulate_Hessian_times_input(DiscretisedDensity<3,elemT>& output,
+    //! Compute the multiplication of the hessian of the prior multiplied by the input.
+    virtual Succeeded accumulate_Hessian_times_input(DiscretisedDensity<3,elemT>& output,
                                                      const DiscretisedDensity<3,elemT>& current_estimate,
                                                      const DiscretisedDensity<3,elemT>& input) const;
 
@@ -226,21 +220,21 @@ class LogcoshPrior:  public
       { return tanh(x)/x; }
     }
 
-    //! The second partial derivatives of the LogCosh Prior
+    //! The Hessian of log(cosh()) is sech^2(x) = (1/cosh(x))^2
     /*!
-     derivative_20 refers to the second derivative w.r.t. x_j only (i.e. diagonal elements of the Hessian)
-     derivative_11 refers to the second derivative w.r.t. x_j and x_k (i.e. off-diagonal elements of the Hessian)
-     * @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
+     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
      */
-    //@{
-    elemT derivative_20(const elemT x_j, const elemT x_k) const;
-    elemT derivative_11(const elemT x_j, const elemT x_k) const;
-    //@}
+    static inline float Hessian(const float d, const float scalar)
+    {
+      const float x = d * scalar;
+      return square((1/ cosh(x)));
+    }
 };
 
 
 END_NAMESPACE_STIR
 
-#endif
+#endif

src/include/stir/recon_buildblock/PLSPrior.h

@@ -121,8 +121,6 @@ class PLSPrior:  public
   /*! \todo set the anatomical image to zero if not defined */
   virtual Succeeded set_up(shared_ptr<const DiscretisedDensity<3,elemT> > const& target_sptr);
 
-  bool is_convex() const;
-
   //! compute the value of the function
   double
     compute_value(const DiscretisedDensity<3,elemT> &current_image_estimate);

src/include/stir/recon_buildblock/QuadraticPrior.h

@@ -103,8 +103,6 @@ class QuadraticPrior:  public
     parabolic_surrogate_curvature_depends_on_argument() const
     { return false; }
 
-  bool is_convex() const;
-
   //! compute the value of the function
   double
     compute_value(const DiscretisedDensity<3,elemT> &current_image_estimate);
@@ -118,20 +116,20 @@ class QuadraticPrior:  public
   void parabolic_surrogate_curvature(DiscretisedDensity<3,elemT>& parabolic_surrogate_curvature, 
                         const DiscretisedDensity<3,elemT> &current_image_estimate);
 
-  virtual void
-  compute_Hessian(DiscretisedDensity<3,elemT>& prior_Hessian_for_single_densel,
-                  const BasicCoordinate<3,int>& coords,
-                  const DiscretisedDensity<3,elemT> &current_image_estimate) const;
+  //! 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);
 
   //! Call accumulate_Hessian_times_input
-  virtual void
+  virtual Succeeded 
     add_multiplication_with_approximate_Hessian(DiscretisedDensity<3,elemT>& output,
                                                 const DiscretisedDensity<3,elemT>& input) const;
 
   //! Compute the multiplication of the hessian of the prior multiplied by the input.
   //! For the quadratic function, the hessian of the prior is 1.
   //! Therefore this will return the weights multiplied by the input.
-  virtual void accumulate_Hessian_times_input(DiscretisedDensity<3,elemT>& output,
+  virtual Succeeded accumulate_Hessian_times_input(DiscretisedDensity<3,elemT>& output,
                                                    const DiscretisedDensity<3,elemT>& current_estimate,
                                                    const DiscretisedDensity<3,elemT>& input) const;
 
@@ -181,19 +179,6 @@ class QuadraticPrior:  public
   virtual bool post_processing();
  private:
   shared_ptr<const DiscretisedDensity<3,elemT> > kappa_ptr;
-
-  //! The second partial derivatives of the Quadratic Prior
-  /*!
-   derivative_20 refers to the second derivative w.r.t. x_j (i.e. diagonal elements of the Hessian)
-   derivative_11 refers to the second derivative w.r.t. x_j and x_k (i.e. off-diagonal elements of the Hessian)
-   * @param x_j is the target voxel.
-   * @param x_k is the voxel in the neighbourhood.
-   * @return the second order partial derivatives of the Quadratic Prior
-   */
-  //@{
-  elemT derivative_20(const elemT x_j, const elemT x_k) const;
-  elemT derivative_11(const elemT x_j, const elemT x_k) const;
-  //@}
 };
 
 

src/include/stir/recon_buildblock/RelativeDifferencePrior.h

@@ -121,19 +121,11 @@ class RelativeDifferencePrior:  public
   void compute_gradient(DiscretisedDensity<3,elemT>& prior_gradient, 
                         const DiscretisedDensity<3,elemT> &current_image_estimate);
 
-  virtual void compute_Hessian(DiscretisedDensity<3,elemT>& prior_Hessian_for_single_densel,
-                                    const BasicCoordinate<3,int>& coords,
-                                    const DiscretisedDensity<3,elemT> &current_image_estimate) const;
 
-  virtual void
+  virtual Succeeded 
     add_multiplication_with_approximate_Hessian(DiscretisedDensity<3,elemT>& output,
                                                 const DiscretisedDensity<3,elemT>& input) const;
 
-    //! Compute the multiplication of the hessian of the prior multiplied by the input.
-  virtual void accumulate_Hessian_times_input(DiscretisedDensity<3,elemT>& output,
-                                                   const DiscretisedDensity<3,elemT>& current_estimate,
-                                                   const DiscretisedDensity<3,elemT>& input) const;
-
   //! get the gamma value used in RDP
   float get_gamma() const;
   //! set the gamma value used in the RDP
@@ -164,8 +156,6 @@ class RelativeDifferencePrior:  public
   //! Has to be called before using this object
   virtual Succeeded set_up(shared_ptr<DiscretisedDensity<3,elemT> > const& target_sptr);
 
-  bool is_convex() const;
-
 protected:
   //! Create variable gamma for Relative Difference Penalty
   float gamma;
@@ -199,22 +189,6 @@ class RelativeDifferencePrior:  public
   virtual bool post_processing();
  private:
   shared_ptr<DiscretisedDensity<3,elemT> > kappa_ptr;
-
-  //! The second partial derivatives of the Relative Difference Prior
-  /*!
-   derivative_20 refers to the second derivative w.r.t. x_j only (i.e. diagonal elements of the Hessian)
-   derivative_11 refers to the second derivative w.r.t. x_j and x_k (i.e. off-diagonal elements of the Hessian)
-   See J. Nuyts, et al., 2002, Equation 7.
-   In the instance x_j, x_k and epsilon equal 0.0, these functions return 0.0 to prevent returning an undefined value
-   due to 0/0 computation. This is a "reasonable" solution to this issue.
-   * @param x_j is the target voxel.
-   * @param x_k is the voxel in the neighbourhood.
-   * @return the second order partial derivatives of the Relative Difference Prior
-   */
-  //@{
-  elemT derivative_20(const elemT x_j, const elemT x_k) const;
-  elemT derivative_11(const elemT x_j, const elemT x_k) const;
-  //@}
 };
 
 

src/recon_buildblock/FilterRootPrior.cxx

@@ -41,12 +41,6 @@ FilterRootPrior(shared_ptr<DataProcessor<DataT> >const& filter_sptr, float penal
   this->penalisation_factor = penalisation_factor_v;
 }
 
-template <typename DataT>
-bool FilterRootPrior<DataT>::
-is_convex() const
-{
-  return false;
-}
 
 template < class T>
 static inline int

src/recon_buildblock/GeneralisedObjectiveFunction.cxx

@@ -263,7 +263,10 @@ add_multiplication_with_approximate_sub_Hessian(TargetT& output,
     {
       // TODO used boost:scoped_ptr
       shared_ptr<TargetT>  prior_output_sptr(output.get_empty_copy());
-      this->prior_sptr->add_multiplication_with_approximate_Hessian(*prior_output_sptr, output);
+      if (this->prior_sptr->add_multiplication_with_approximate_Hessian(*prior_output_sptr, output) ==
+	  Succeeded::no)
+	return Succeeded::no;
+
 
        // (*prior_output_sptr)/= num_subsets;
        // output -= *prior_output_sptr;
@@ -380,7 +383,9 @@ accumulate_sub_Hessian_times_input(TargetT& output,
   {
     // TODO used boost:scoped_ptr
     shared_ptr<TargetT>  prior_output_sptr(output.get_empty_copy());
-    this->prior_sptr->accumulate_Hessian_times_input(*prior_output_sptr, current_image_estimate, output);
+    if (this->prior_sptr->accumulate_Hessian_times_input(*prior_output_sptr, current_image_estimate, output) ==
+        Succeeded::no)
+      return Succeeded::no;
 
     typename TargetT::const_full_iterator prior_output_iter = prior_output_sptr->begin_all_const();
     const typename TargetT::const_full_iterator end_prior_output_iter = prior_output_sptr->end_all_const();

src/recon_buildblock/GeneralisedPrior.cxx

@@ -52,37 +52,26 @@ set_up(shared_ptr<const TargetT> const&)
 }
 
 template <typename TargetT>
-void
-GeneralisedPrior<TargetT>::
-compute_Hessian(TargetT& output,
-                const BasicCoordinate<3,int>& coords,
-                const TargetT& current_image_estimate) const
-{
-  if (this->is_convex())
-    error("GeneralisedPrior:\n  compute_Hessian implementation is not overloaded by your convex prior.");
-  else
-    error("GeneralisedPrior:\n  compute_Hessian is not implemented for this (non-convex) prior.");
-}
-
-template <typename TargetT>
-void
+Succeeded 
 GeneralisedPrior<TargetT>::
 add_multiplication_with_approximate_Hessian(TargetT& output,
               const TargetT& input) const
 {
   error("GeneralisedPrior:\n"
   "add_multiplication_with_approximate_Hessian implementation is not overloaded by your prior.");
+  return Succeeded::no;
 }
 
 template <typename TargetT>
-void
+Succeeded
 GeneralisedPrior<TargetT>::
 accumulate_Hessian_times_input(TargetT& output,
         const TargetT& current_estimate,
         const TargetT& input) const
 {
   error("GeneralisedPrior:\n"
         "accumulate_Hessian_times_input implementation is not overloaded by your prior.");
+  return Succeeded::no;
 }
 
 template <typename TargetT>