Skip to content

Commit ffd3ff8

Browse files
author
robbietuk
authored
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
1 parent 9d93017 commit ffd3ff8

17 files changed

Lines changed: 301 additions & 126 deletions

‎documentation/release_5.0.htm‎

Lines changed: 6 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -143,6 +143,12 @@ <h3>Changed functionality</h3>
143143
<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.
144144
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).
145145
Futhermore, added functionality to exclude true and unscattered event types from list mode processing. </li>
146+
<li>
147+
Poisson log-likelihood gradient methods now use <tt>actual_compute_subset_gradient_without_penalty</tt>
148+
to properly handle the objective function gradient computation at the distributable level.
149+
Removes the subtraction of the subset sensitivity, which lead to a inconsistency between log-likelihood function
150+
and gradient when <tt>use_subset_sensitivities</tt> was false.
151+
</li>
146152
</ul>
147153

148154

‎recon_test_pack/OSSPS_test_PM_QP_subsens1.par‎

Lines changed: 73 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -0,0 +1,73 @@
1+
OSSPSParameters :=
2+
; sample file for OSSPS
3+
; parameters used here are for illustrative purposes only
4+
; i.e. they are not recommended values
5+
6+
objective function type:= PoissonLogLikelihoodWithLinearModelForMeanAndProjData
7+
PoissonLogLikelihoodWithLinearModelForMeanAndProjData Parameters:=
8+
9+
input file := Utahscat600k_ca_seg4.hs
10+
zero end planes of segment 0:= 1
11+
; if disabled, defaults to maximum segment number in the file
12+
maximum absolute segment number to process := 3
13+
14+
; change to STIR 2.x default for compatibility
15+
use subset sensitivities:=1
16+
sensitivity filename:= RPTsens_seg3_PM.hv
17+
18+
projector pair type := Matrix
19+
Projector Pair Using Matrix Parameters :=
20+
Matrix type := Ray Tracing
21+
Ray tracing matrix parameters :=
22+
number of rays in tangential direction to trace for each bin := 2
23+
; restrict to cylindrical fov := 0
24+
End Ray tracing matrix parameters :=
25+
End Projector Pair Using Matrix Parameters :=
26+
27+
; additive sinogram:=my_fake_randoms.hs
28+
prior type := quadratic
29+
Quadratic Prior Parameters:=
30+
penalisation factor := 0.5
31+
; next defaults to 0, set to 1 for 2D inverse Euclidean weights, 0 for 3D
32+
only 2D:= 0
33+
END Quadratic Prior Parameters:=
34+
35+
End PoissonLogLikelihoodWithLinearModelForMeanAndProjData Parameters:=
36+
37+
output filename prefix := my_test_image_OSSPS_PM_QP_subsens
38+
; iteration scheme
39+
40+
number of subsets:= 4
41+
;start at subset:= 0
42+
;start at subiteration number := 1
43+
number of subiterations:= 8
44+
Save estimates at subiteration intervals:= 8
45+
;write update image := 0
46+
;report objective function values interval := 1
47+
48+
; if next is disabled, defaults to image full of 1s (but that's not good for OSSPS)
49+
; in particular, make sure it has the correct scale
50+
initial estimate:= test_image_PM_QP_6.hv
51+
enforce initial positivity condition := 1
52+
53+
; here start OSSPS specific values
54+
55+
; values to use for the 'precomputed denominator'
56+
; specify either procomputed denomiator or normalisation type
57+
; use the following if you have it already (e.g. from previous run)
58+
; note: setting the value to 1 will use an images full of ones.
59+
; precomputed denominator := my_precomputed_denominator.hv
60+
61+
; specify relaxation scheme
62+
; lambda = relaxation_parameter/ (1+relaxation_gamma*(subiteration_num/num_subsets)
63+
relaxation parameter := 2
64+
relaxation gamma:=.1
65+
66+
67+
END :=
68+
69+
70+
71+
72+
73+

‎recon_test_pack/run_tests.sh‎

Lines changed: 15 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -245,6 +245,21 @@ echo There were problems here!;
245245
ThereWereErrors=1;
246246
fi
247247

248+
${MPIRUN} ${INSTALL_DIR}OSSPS OSSPS_test_PM_QP_subsens1.par 1> OSSPS_PM_QP.log 2> OSSPS_PM_QP_stderr.log
249+
250+
echo '---- Comparing output of OSSPS subiter 8 using subset sensitivity (should be identical up to tolerance)'
251+
echo Running ${INSTALL_DIR}compare_image
252+
# relax test for the outer-rim voxels as these turn out to be more unstable than the internal ones
253+
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
254+
${INSTALL_DIR}compare_image -r 1 test_image_OSSPS_PM_QP_8.hv my_test_image_OSSPS_PM_QP_subsens_8.hv
255+
then
256+
echo ---- This test seems to be ok !;
257+
else
258+
echo There were problems here!;
259+
ThereWereErrors=1;
260+
fi
261+
262+
248263
echo
249264
echo ------------- tests on stir_math and correct_projdata ---------
250265
echo "first make up some randoms (just a projdata full of 1)"

‎recon_test_pack/test_image_OSSPS_PM_QP_8.v‎

0 Bytes
Binary file not shown.

‎src/include/stir/recon_buildblock/PoissonLogLikelihoodWithLinearKineticModelAndDynamicProjectionData.h‎

Lines changed: 7 additions & 4 deletions
Original file line numberDiff line numberDiff line change
@@ -70,10 +70,6 @@ public RegisteredParsingObject<PoissonLogLikelihoodWithLinearKineticModelAndDyn
7070
*/
7171
virtual TargetT *
7272
construct_target_ptr() const;
73-
virtual void
74-
compute_sub_gradient_without_penalty_plus_sensitivity(TargetT& gradient,
75-
const TargetT &current_estimate,
76-
const int subset_num);
7773

7874
virtual std::unique_ptr<ExamInfo>
7975
get_exam_info_uptr_for_target() const;
@@ -170,6 +166,13 @@ public RegisteredParsingObject<PoissonLogLikelihoodWithLinearKineticModelAndDyn
170166

171167
bool actual_subsets_are_approximately_balanced(std::string& warning_message) const;
172168

169+
170+
virtual void
171+
actual_compute_subset_gradient_without_penalty(TargetT& gradient,
172+
const TargetT &current_estimate,
173+
const int subset_num,
174+
const bool add_sensitivity);
175+
173176
//! Sets defaults for parsing
174177
/*! Resets \c sensitivity_filename and \c sensitivity_sptr and
175178
\c recompute_sensitivity to \c false.

‎src/include/stir/recon_buildblock/PoissonLogLikelihoodWithLinearKineticModelAndDynamicProjectionData.txx‎

Lines changed: 8 additions & 6 deletions
Original file line numberDiff line numberDiff line change
@@ -426,9 +426,10 @@ set_normalisation_sptr(const shared_ptr<BinNormalisation>& arg)
426426
template<typename TargetT>
427427
void
428428
PoissonLogLikelihoodWithLinearKineticModelAndDynamicProjectionData<TargetT>::
429-
compute_sub_gradient_without_penalty_plus_sensitivity(TargetT& gradient,
430-
const TargetT &current_estimate,
431-
const int subset_num)
429+
actual_compute_subset_gradient_without_penalty(TargetT& gradient,
430+
const TargetT &current_estimate,
431+
const int subset_num,
432+
const bool add_sensitivity)
432433
{
433434
if (subset_num<0 || subset_num>=this->get_num_subsets())
434435
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,
452453

453454

454455
this->_single_frame_obj_funcs[frame_num].
455-
compute_sub_gradient_without_penalty_plus_sensitivity(dyn_gradient[frame_num],
456-
dyn_image_estimate[frame_num],
457-
subset_num);
456+
actual_compute_subset_gradient_without_penalty(dyn_gradient[frame_num],
457+
dyn_image_estimate[frame_num],
458+
subset_num,
459+
add_sensitivity);
458460
}
459461

460462
this->_patlak_plot_sptr->multiply_dynamic_image_with_model_gradient(gradient,

‎src/include/stir/recon_buildblock/PoissonLogLikelihoodWithLinearModelForMean.h‎

Lines changed: 54 additions & 32 deletions
Original file line numberDiff line numberDiff line change
@@ -66,14 +66,13 @@ START_NAMESPACE_STIR
6666
the <i>sensitivity</i> because (if \f$r=0\f$) it is the total
6767
probability of detecting a count (in any bin) originating from \f$v\f$.
6868
69-
This class computes the gradient as a sum of these two terms. The
70-
sensitivity has to be computed by the virtual function
71-
\c add_subset_sensitivity(). The sum is computed by
72-
\c compute_sub_gradient_without_penalty_plus_sensitivity().
73-
74-
The reason for this is that the sensitivity is data-independent, and
75-
can be computed only once. See also
76-
PoissonLogLikelihoodWithLinearModelForMeanAndListModeData.
69+
This class computes the gradient directly, via \c compute_sub_gradient_without_penalty().
70+
This method is utilised by the \c OSSPS algorithm in STIR.
71+
However, an additional method (\c compute_sub_gradient_without_penalty_plus_sensitivity())
72+
is provided that computes the sum of the subset gradient (without penalty) and the sensitivity.
73+
This method is utilised by the \c OSMAPOSL algorithm.
74+
75+
See also \c PoissonLogLikelihoodWithLinearModelForMeanAndListModeData.
7776
7877
\par Relation with Kullback-Leibler distance
7978
@@ -120,38 +119,40 @@ public GeneralisedObjectiveFunction<TargetT>
120119

121120
//PoissonLogLikelihoodWithLinearModelForMean();
122121

123-
//! Implementation in terms of compute_sub_gradient_without_penalty_plus_sensitivity()
124-
/*! \warning If separate subsensitivities are not used, we just subtract the total
125-
sensitivity divided by the number of subsets.
126-
This is fine for some algorithms as the sum over all the subsets is
127-
equal to gradient of the objective function (without prior).
128-
Other algorithms do not behave very stable under this approximation
129-
however. So, currently setup() will return an error if
130-
<code>!subsets_are_approximately_balanced()</code> and subset sensitivities
131-
are not used.
132-
133-
\see get_use_subset_sensitivities()
134-
*/
122+
//! Compute the subset gradient of the (unregularised) objective function
123+
/*!
124+
Implementation in terms of actual_compute_sub_gradient_without_penalty()
125+
This function is used by OSSPS may be used by other gradient ascent/descent algorithms
126+
127+
This computes
128+
\f[
129+
{\partial L \over \partial \lambda_v} =
130+
\sum_b P_{bv} ({y_b \over Y_b} - 1)
131+
\f]
132+
(see the class general documentation).
133+
The sum will however be restricted to a subset.
134+
*/
135135
virtual void
136136
compute_sub_gradient_without_penalty(TargetT& gradient,
137137
const TargetT &current_estimate,
138138
const int subset_num);
139139

140-
//! This should compute the gradient of the (unregularised) objective function plus the (sub)sensitivity
141-
/*!
142-
This function is used for instance by OSMAPOSL.
143-
144-
This computes
145-
\f[ {\partial L \over \partial \lambda_v} + P_v =
146-
\sum_b P_{bv} {y_b \over Y_b}
147-
\f]
148-
(see the class general documentation).
149-
The sum will however be restricted to a subset.
150-
*/
140+
//! This should compute the subset gradient of the (unregularised) objective function plus the subset sensitivity
141+
/*!
142+
Implementation in terms of actual_compute_sub_gradient_without_penalty().
143+
This function is used for instance by OSMAPOSL.
144+
145+
This computes
146+
\f[ {\partial L \over \partial \lambda_v} + P_v =
147+
\sum_b P_{bv} {y_b \over Y_b}
148+
\f]
149+
(see the class general documentation).
150+
The sum will however be restricted to a subset.
151+
*/
151152
virtual void
152153
compute_sub_gradient_without_penalty_plus_sensitivity(TargetT& gradient,
153154
const TargetT &current_estimate,
154-
const int subset_num) =0;
155+
const int subset_num);
155156

156157
//! set-up sensitivity etc if possible
157158
/*! If \c recompute_sensitivity is \c false, we will try to
@@ -264,6 +265,27 @@ public GeneralisedObjectiveFunction<TargetT>
264265
*/
265266
void compute_sensitivities();
266267

268+
//! computes the subset gradient of the objective function without the penalty (optional: add subset sensitivity)
269+
/*!
270+
If \c add_sensitivity is \c true, this computes
271+
\f[ {\partial L \over \partial \lambda_v} + P_v =
272+
\sum_b P_{bv} {y_b \over Y_b}
273+
\f]
274+
(see the class general documentation).
275+
The sum will however be restricted to a subset.
276+
277+
However, if \c add_sensitivity is \c false, this function will instead compute only the gradient
278+
\f[
279+
{\partial L \over \partial \lambda_v} =
280+
\sum_b P_{bv} ({y_b \over Y_b} - 1)
281+
\f]
282+
*/
283+
virtual void
284+
actual_compute_subset_gradient_without_penalty(TargetT& gradient,
285+
const TargetT &current_estimate,
286+
const int subset_num,
287+
const bool add_sensitivity) = 0;
288+
267289
//! Sets defaults for parsing
268290
/*! Resets \c sensitivity_filename, \c subset_sensitivity_filenames to empty,
269291
\c recompute_sensitivity to \c false, and \c use_subset_sensitivities to false.

‎src/include/stir/recon_buildblock/PoissonLogLikelihoodWithLinearModelForMeanAndGatedProjDataWithMotion.h‎

Lines changed: 5 additions & 4 deletions
Original file line numberDiff line numberDiff line change
@@ -79,10 +79,11 @@ public RegisteredParsingObject<PoissonLogLikelihoodWithLinearModelForMeanAndGat
7979
virtual TargetT *
8080
construct_target_ptr() const;
8181

82-
virtual void
83-
compute_sub_gradient_without_penalty_plus_sensitivity(TargetT& gradient,
84-
const TargetT &current_estimate,
85-
const int subset_num);
82+
virtual void
83+
actual_compute_subset_gradient_without_penalty(TargetT& gradient,
84+
const TargetT &current_estimate,
85+
const int subset_num,
86+
const bool add_sensitivity);
8687

8788
virtual double
8889
actual_compute_objective_function_without_penalty(const TargetT& current_estimate,

‎src/include/stir/recon_buildblock/PoissonLogLikelihoodWithLinearModelForMeanAndGatedProjDataWithMotion.txx‎

Lines changed: 8 additions & 6 deletions
Original file line numberDiff line numberDiff line change
@@ -415,9 +415,10 @@ set_up_before_sensitivity(shared_ptr<const TargetT > const& target_sptr)
415415
template<typename TargetT>
416416
void
417417
PoissonLogLikelihoodWithLinearModelForMeanAndGatedProjDataWithMotion<TargetT>::
418-
compute_sub_gradient_without_penalty_plus_sensitivity(TargetT& gradient,
419-
const TargetT &current_estimate,
420-
const int subset_num)
418+
actual_compute_subset_gradient_without_penalty(TargetT& gradient,
419+
const TargetT &current_estimate,
420+
const int subset_num,
421+
const bool add_sensitivity)
421422
{
422423
assert(subset_num>=0);
423424
assert(subset_num<this->num_subsets);
@@ -435,9 +436,10 @@ compute_sub_gradient_without_penalty_plus_sensitivity(TargetT& gradient,
435436
gated_gradient[gate_num].end_all(),
436437
0.F);
437438
this->_single_gate_obj_funcs[gate_num].
438-
compute_sub_gradient_without_penalty_plus_sensitivity(gated_gradient[gate_num],
439-
gated_image_estimate[gate_num],
440-
subset_num);
439+
actual_compute_subset_gradient_without_penalty(gated_gradient[gate_num],
440+
gated_image_estimate[gate_num],
441+
subset_num,
442+
add_sensitivity);
441443
}
442444
// if(this->_motion_correction_type==-1)
443445
this->_reverse_motion_vectors.warp_image(gradient,gated_gradient) ;

‎src/include/stir/recon_buildblock/PoissonLogLikelihoodWithLinearModelForMeanAndListModeData.h‎

Lines changed: 0 additions & 4 deletions
Original file line numberDiff line numberDiff line change
@@ -78,10 +78,6 @@ public PoissonLogLikelihoodWithLinearModelForMean<TargetT>
7878
virtual Succeeded
7979
set_up(shared_ptr <TargetT > const& target_sptr);
8080

81-
virtual
82-
void compute_sub_gradient_without_penalty_plus_sensitivity(TargetT& gradient,
83-
const TargetT &current_estimate,
84-
const int subset_num)=0;
8581
//! time frame definitions
8682
/*! \todo This is currently used to be able to compute the gradient for
8783
one time frame. However, it probably does not belong here.

0 commit comments

Comments
 (0)