From 2992bc589cdfab7f44f7b92f2978ad7b96a5ac89 Mon Sep 17 00:00:00 2001
From: Kris Thielemans <k.thielemans@ucl.ac.uk>
Date: Mon, 22 Jul 2024 22:58:39 +0100
Subject: [PATCH] Hessian times input: throw error if result would be incorrect

The Hessian calculations for the Poisson log-lik for projdata
use divide_and_truncate, which sets negative
numerators to zero. This is incorrect when the forward projection
of the "input" contains negatives.
We now check this and throw an error if it occurs. A proper fix
will have to be for later.

See https://github.com/UCL/STIR/issues/1461
---
 documentation/release_6.2.htm                 |  7 +++++
 ...ihoodWithLinearModelForMeanAndProjData.cxx | 26 +++++++++++++++++++
 2 files changed, 33 insertions(+)

diff --git a/documentation/release_6.2.htm b/documentation/release_6.2.htm
index 6bc503297..7b544a000 100644
--- a/documentation/release_6.2.htm
+++ b/documentation/release_6.2.htm
@@ -106,6 +106,13 @@ <h3>Bug fixes</h3>
     handled. Also, downsampling of <code>BlocksOnCylindrical</code> scanners in scatter simulation was inaccurate.<br>
     <a href=https://github.com/UCL/STIR/pull/1466>PR #1466</a>
   </li>
+  <li>
+    The "Hessian times input" calculations of the Poisson log-likelihood for projection data
+    were incorrect when the forward projection of the "input" contains negatives.
+    We now detect this and throw an error if it occurs. A proper fix
+    will have to be for later.<br>
+    See <a href="https://github.com/UCL/STIR/issues/1461">Issue #1461</a>
+  </li>
 </ul>
 
 <h3>Build system</h3>
diff --git a/src/recon_buildblock/PoissonLogLikelihoodWithLinearModelForMeanAndProjData.cxx b/src/recon_buildblock/PoissonLogLikelihoodWithLinearModelForMeanAndProjData.cxx
index 9dd07530d..e6c6cdd5d 100644
--- a/src/recon_buildblock/PoissonLogLikelihoodWithLinearModelForMeanAndProjData.cxx
+++ b/src/recon_buildblock/PoissonLogLikelihoodWithLinearModelForMeanAndProjData.cxx
@@ -977,12 +977,15 @@ PoissonLogLikelihoodWithLinearModelForMeanAndProjData<
                                               this->get_num_subsets());
 
   info("Forward projecting input image.", 2);
+  volatile bool any_negatives = false;
 #ifdef STIR_OPENMP
 #  pragma omp parallel for schedule(dynamic)
 #endif
   // note: older versions of openmp need an int as loop
   for (int i = 0; i < static_cast<int>(vg_idx_to_process.size()); ++i)
     {
+      if (any_negatives)
+        continue; // early exit as we'll throw error outside of the parallel for
       const auto viewgram_idx = vg_idx_to_process[i];
       {
 #ifdef STIR_OPENMP
@@ -1011,6 +1014,11 @@ PoissonLogLikelihoodWithLinearModelForMeanAndProjData<
       {
         tmp_viewgrams = this->get_proj_data().get_empty_related_viewgrams(viewgram_idx, symmetries_sptr);
         this->get_projector_pair().get_forward_projector_sptr()->forward_project(tmp_viewgrams);
+        if (tmp_viewgrams.find_min() < 0)
+          {
+            any_negatives = true;
+            continue; // throw error outside of parallel for
+          }
       }
 
       // now divide by the data term
@@ -1024,6 +1032,11 @@ PoissonLogLikelihoodWithLinearModelForMeanAndProjData<
 
     } // end of loop over view/segments
 
+  if (any_negatives)
+    error("PoissonLL add_multiplication_with_approximate_sub_Hessian: forward projection of input contains negatives. The "
+          "result would be incorrect, so we abort.\n"
+          "See https://github.com/UCL/STIR/issues/1461");
+
   shared_ptr<TargetT> tmp(output.get_empty_copy());
   this->get_projector_pair().get_back_projector_sptr()->get_output(*tmp);
   // output += tmp;
@@ -1098,11 +1111,15 @@ PoissonLogLikelihoodWithLinearModelForMeanAndProjData<TargetT>::actual_accumulat
 
   // Forward project input image
   info("Forward projecting input image.", 2);
+  volatile bool any_negatives = false;
 #ifdef STIR_OPENMP
 #  pragma omp parallel for schedule(dynamic)
 #endif
   for (int i = 0; i < static_cast<int>(vg_idx_to_process.size()); ++i)
     { // Loop over each of the viewgrams in input_viewgrams_vec, forward projecting input into them
+      if (any_negatives)
+        continue; // early exit as we'll throw error outside of the parallel for
+
       const auto viewgram_idx = vg_idx_to_process[i];
       {
 #ifdef STIR_OPENMP
@@ -1118,7 +1135,16 @@ PoissonLogLikelihoodWithLinearModelForMeanAndProjData<TargetT>::actual_accumulat
       }
       input_viewgrams_vec[i] = this->get_proj_data().get_empty_related_viewgrams(viewgram_idx, symmetries_sptr);
       this->get_projector_pair().get_forward_projector_sptr()->forward_project(input_viewgrams_vec[i]);
+      if (input_viewgrams_vec[i].find_min() < 0)
+        {
+          any_negatives = true;
+          continue; // throw error outside of parallel for
+        }
     }
+  if (any_negatives)
+    error("PoissonLL accumulate_sub_Hessian_times_input: forward projection of input contains negatives. The "
+          "result would be incorrect, so we abort.\n"
+          "See https://github.com/UCL/STIR/issues/1461");
 
   info("Forward projecting current image estimate and back projecting to output.", 2);
   this->get_projector_pair().get_forward_projector_sptr()->set_input(current_image_estimate);