STIR 6.4.0
stir::PoissonLogLikelihoodWithLinearModelForMean< TargetT > Class Template Referenceabstract

a base class for LogLikelihood of independent Poisson variables where the mean values are linear combinations of the parameters. More...

#include "stir/recon_buildblock/PoissonLogLikelihoodWithLinearModelForMean.h"

Inheritance diagram for stir::PoissonLogLikelihoodWithLinearModelForMean< TargetT >:

Public Member Functions

void compute_sub_gradient_without_penalty (TargetT &gradient, const TargetT &current_estimate, const int subset_num) override
 Compute the subset gradient of the (unregularised) objective function.
 
virtual void compute_sub_gradient_without_penalty_plus_sensitivity (TargetT &gradient, const TargetT &current_estimate, const int subset_num)
 This should compute the subset gradient of the (unregularised) objective function plus the subset sensitivity.
 
Succeeded set_up (shared_ptr< TargetT > const &target_sptr) override
 set-up sensitivity etc if possible
 
const TargetT & get_sensitivity () const
 Get a const reference to the total sensitivity.
 
const TargetT & get_subset_sensitivity (const int subset_num) const
 Get a const reference to the sensitivity for a subset.
 
virtual void add_subset_sensitivity (TargetT &sensitivity, const int subset_num) const =0
 Add subset sensitivity to existing data.
 
bool get_use_subset_sensitivities () const
 find out if subset_sensitivities are used
 
bool get_recompute_sensitivity () const
 find current value of recompute_sensitivity
 
std::string get_sensitivity_filename () const
 get filename to read (or write) the total sensitivity
 
std::string get_subsensitivity_filenames () const
 get filename pattern to read (or write) the subset sensitivities
 
std::string get_subsensitivity_filename (const int subset_num) const
 Return the filename for a particular subset.
 
- Public Member Functions inherited from stir::GeneralisedObjectiveFunction< TargetT >
virtual TargetT * construct_target_ptr () const =0
 Creates a suitable target as determined by the parameters.
 
virtual void compute_sub_gradient (TargetT &gradient, const TargetT &current_estimate, const int subset_num)
 Compute the subset-gradient of the objective function at current_estimate.
 
virtual void compute_gradient (TargetT &gradient, const TargetT &current_estimate)
 Compute the gradient of the objective function at the current_estimate.
 
virtual void compute_gradient_without_penalty (TargetT &gradient, const TargetT &current_estimate)
 Compute the gradient of the unregularised objective function at the current_estimate.
 
virtual double compute_objective_function_without_penalty (const TargetT &current_estimate, const int subset_num)
 Compute the value of the unregularised sub-objective function at the current_estimate.
 
virtual double compute_objective_function_without_penalty (const TargetT &current_estimate)
 Compute the value of the unregularised objective function at the current_estimate.
 
double compute_penalty (const TargetT &current_estimate, const int subset_num)
 Compute the value of the sub-penalty at the current_estimate.
 
double compute_penalty (const TargetT &current_estimate)
 Compute the value of the penalty at the current_estimate.
 
double compute_objective_function (const TargetT &current_estimate, const int subset_num)
 Compute the value of the sub-objective function at the current_estimate.
 
double compute_objective_function (const TargetT &current_estimate)
 Compute the value of the objective function at the current_estimate.
 
double compute_value (const TargetT &current_estimate)
 Alias for compute_objective_function(const TargetT&)
 
Succeeded add_multiplication_with_approximate_sub_Hessian_without_penalty (TargetT &output, const TargetT &input, const int subset_num) const
 Functions that multiply the approximate (sub)Hessian with a \'vector\'.
 
Succeeded add_multiplication_with_approximate_sub_Hessian (TargetT &output, const TargetT &input, const int subset_num) const
 
Succeeded add_multiplication_with_approximate_Hessian_without_penalty (TargetT &output, const TargetT &input) const
 
Succeeded add_multiplication_with_approximate_Hessian (TargetT &output, const TargetT &input) const
 
Succeeded accumulate_Hessian_times_input (TargetT &output, const TargetT &current_image_estimate, const TargetT &input) const
 Functions that multiply the True (sub)Hessian with a \'vector\'.
 
Succeeded accumulate_Hessian_times_input_without_penalty (TargetT &output, const TargetT &current_image_estimate, const TargetT &input) const
 
Succeeded accumulate_sub_Hessian_times_input (TargetT &output, const TargetT &current_image_estimate, const TargetT &input, const int subset_num) const
 
Succeeded accumulate_sub_Hessian_times_input_without_penalty (TargetT &output, const TargetT &current_image_estimate, const TargetT &input, const int subset_num) const
 
std::string get_objective_function_values_report (const TargetT &current_estimate)
 Construct a string with info on the value of objective function with and without penalty.
 
int get_num_subsets () const
 Return the number of subsets in-use.
 
virtual std::unique_ptr< ExamInfoget_exam_info_uptr_for_target () const
 
virtual int set_num_subsets (const int num_subsets)=0
 Attempts to change the number of subsets.
 
bool subsets_are_approximately_balanced () const
 Checks of the current subset scheme is approximately balanced.
 
bool subsets_are_approximately_balanced (std::string &warning_message) const
 Checks of the current subset scheme is approximately balanced and constructs a warning message.
 
bool prior_is_zero () const
 check if the prior is set (or the penalisation factor is 0)
 
GeneralisedPrior< TargetT > *const get_prior_ptr () const
 Read-only access to the prior.
 
shared_ptr< GeneralisedPrior< TargetT > > get_prior_sptr ()
 
void set_prior_sptr (const shared_ptr< GeneralisedPrior< TargetT > > &)
 Change the prior.
 
virtual void set_input_data (const shared_ptr< ExamData > &)=0
 set_input_data
 
virtual const ExamDataget_input_data () const =0
 get input data
 
virtual void set_additive_proj_data_sptr (const shared_ptr< ExamData > &)=0
 set_additive_proj_data_sptr
 
virtual void set_normalisation_sptr (const shared_ptr< BinNormalisation > &)=0
 set_normalisation_sptr
 
- Public Member Functions inherited from stir::RegisteredObjectBase
virtual std::string get_registered_name () const =0
 Returns the name of the type of the object.
 
- Public Member Functions inherited from stir::ParsingObject
 ParsingObject (const ParsingObject &)
 
ParsingObjectoperator= (const ParsingObject &)
 
bool parse (std::istream &f)
 
bool parse (const char *const filename)
 
void ask_parameters ()
 
virtual std::string parameter_info ()
 

Functions to set parameters

This can be used as alternative to the parsing mechanism.

Warning
After using any of these, you have to call set_up().
Be careful with setting shared pointers. If you modify the objects in one place, all objects that use the shared pointer will be affected.
void set_recompute_sensitivity (const bool)
 
void set_subset_sensitivity_sptr (const shared_ptr< TargetT > &, const int subset_num)
 
void set_use_subset_sensitivities (const bool)
 See get_use_subset_sensitivities()
 
void set_sensitivity_filename (const std::string &)
 set filename to read (or write) the total sensitivity
 
void set_subsensitivity_filenames (const std::string &)
 set filename pattern to read (or write) the subset sensitivities
 
void fill_nonidentifiable_target_parameters (TargetT &target, const float value) const override
 
virtual Succeeded set_up_before_sensitivity (shared_ptr< const TargetT > const &target_sptr)=0
 set-up specifics for the derived class
 
void compute_sensitivities ()
 compute subset and total sensitivity
 
virtual void actual_compute_subset_gradient_without_penalty (TargetT &gradient, const TargetT &current_estimate, const int subset_num, const bool add_sensitivity)=0
 computes the subset gradient of the objective function without the penalty (optional: add subset sensitivity)
 
void set_defaults () override
 Sets defaults for parsing.
 
void initialise_keymap () override
 sets parsing keys
 
bool post_processing () override
 This will be called at the end of the parsing.
 

Additional Inherited Members

- Static Public Member Functions inherited from stir::RegisteredObject< GeneralisedObjectiveFunction< TargetT > >
static GeneralisedObjectiveFunction< TargetT > * read_registered_object (std::istream *in, const std::string &registered_name)
 Construct a new object (of a type derived from Root, its actual type determined by the registered_name parameter) by parsing the istream.
 
static GeneralisedObjectiveFunction< TargetT > * ask_type_and_parameters ()
 ask the user for the type, and then calls read_registered_object(0, type)
 
static void list_registered_names (std::ostream &stream)
 List all possible registered names to the stream.
 
- Protected Types inherited from stir::RegisteredObject< GeneralisedObjectiveFunction< TargetT > >
typedef GeneralisedObjectiveFunction< TargetT > *(* RootFactory) (std::istream *)
 The type of a root factory is a function, taking an istream* as argument, and returning a Root*.
 
typedef FactoryRegistry< std::string, RootFactory, interfile_lessRegistryType
 The type of the registry.
 
virtual bool actual_subsets_are_approximately_balanced (std::string &warning_message) const =0
 Implementation of function that checks subset balancing.
 
virtual double actual_compute_objective_function_without_penalty (const TargetT &current_estimate, const int subset_num)=0
 Implementation of function that computes the objective function for the current subset.
 
virtual Succeeded actual_add_multiplication_with_approximate_sub_Hessian_without_penalty (TargetT &output, const TargetT &input, const int subset_num) const
 Implementation of the function that multiplies the approximate sub-Hessian with a vector.
 
virtual Succeeded actual_accumulate_sub_Hessian_times_input_without_penalty (TargetT &output, const TargetT &current_image_estimate, const TargetT &input, const int subset_num) const
 Implementation of the function computes the sub-Hessian and multiplies by a vector.
 
virtual void set_key_values ()
 This will be called before parsing or parameter_info is called.
 
- Static Protected Member Functions inherited from stir::RegisteredObject< GeneralisedObjectiveFunction< TargetT > >
static RegistryTyperegistry ()
 Static function returning the registry.
 
int num_subsets
 
bool already_set_up
 
shared_ptr< GeneralisedPrior< TargetT > > prior_sptr
 
KeyParser parser
 

Detailed Description

template<typename TargetT>
class stir::PoissonLogLikelihoodWithLinearModelForMean< TargetT >

a base class for LogLikelihood of independent Poisson variables where the mean values are linear combinations of the parameters.

Suppose the measured data $y$ are statistically independent and Poisson distributed with mean $\bar y$. The log likelihood is then given by

\[ L(y\;\bar y) = \sum_b y_b log \bar y_b - \mathrm{log} y_b! - \bar y_b \]

If the means are modeled by some functions of some parameters $Y(\lambda)$, the gradient w.r.t. those parameters follows immediately as

\[ {\partial L \over \partial \lambda_v} =
    \sum_b \left( {y_b \over Y_b} - 1 \right)
           {\partial Y_b \over \partial \lambda_v }
\]

In this class, it is assumed that the means are linear functions of the parameters, which using matrices and vectors can be written (in a notation familiar to PET people) as

\[ Y = P \lambda + r \]

If the background term $r$ is 0, $P_{bv}$ gives the conditional probability of detecting a count in bin $b$ if the parameter $\lambda_v$ is 1 and all others 0.

In this case, the gradient is given by

\[ {\partial L \over \partial \lambda_v} =
    \sum_b P_{bv} {y_b \over Y_b} - P_v
\]

where

\[P_v = \sum_b P_{bv}
\]

where the sum is over all possible bins (not only those where any counts were detected). We call the vector given by $P_v$ the sensitivity because (if $r=0$) it is the total probability of detecting a count (in any bin) originating from $v$.

This class computes the gradient directly, via compute_sub_gradient_without_penalty(). This method is utilised by the OSSPS algorithm in STIR. However, an additional method (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 OSMAPOSL algorithm.

See also PoissonLogLikelihoodWithLinearModelForMeanAndListModeData.

Relation with Kullback-Leibler distance

Note that up to terms independent of $\bar y$, the log likelihood is equal to minus the Kullback-Leibler generalised distance

\[ \mathrm{KL}(y,\bar y) = \sum_b y_b \mathrm{log}(y_b/\bar y_b) + \bar y_b - y_b \]

which has the nice property that $\mathrm{KL}(y,\bar y)=0$ implies that $y=\bar y$.

Parameters for parsing
Defaults are indicated below
; specifies if we keep separate sensitivity images (which is more accurate and is
; recommended) or if we assume the subsets are exactly balanced
; and hence compute the subset-senstivity as sensitivity/num_subsets (this uses less memory).
use_subset_sensitivities := 1
; for recomputing sensitivity, even if a filename is specified
recompute sensitivity:= 0
; filename for reading the sensitivity, or writing if it is recomputed
; if use_subset_sensitivities=0
sensitivity filename:=
; pattern for filename for reading the subset sensitivities, or writing if recomputed
; if use_subset_sensitivities=1
; e.g. subsens_%d.hv
; fmt::format is used with the pattern
subset sensitivity filenames:=
Terminology
We currently use sub_gradient for the gradient of the likelihood of the subset (not the mathematical subgradient).

Member Function Documentation

◆ compute_sub_gradient_without_penalty()

template<typename TargetT>
void stir::PoissonLogLikelihoodWithLinearModelForMean< TargetT >::compute_sub_gradient_without_penalty ( TargetT & gradient,
const TargetT & current_estimate,
const int subset_num )
overridevirtual

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

\[ {\partial L \over \partial \lambda_v} =
   \sum_b P_{bv} ({y_b \over Y_b} - 1)
\]

(see the class general documentation). The sum will however be restricted to a subset.

Implements stir::GeneralisedObjectiveFunction< TargetT >.

References actual_compute_subset_gradient_without_penalty(), compute_sub_gradient_without_penalty(), and stir::error().

Referenced by compute_sub_gradient_without_penalty().

◆ compute_sub_gradient_without_penalty_plus_sensitivity()

template<typename TargetT>
void stir::PoissonLogLikelihoodWithLinearModelForMean< TargetT >::compute_sub_gradient_without_penalty_plus_sensitivity ( TargetT & gradient,
const TargetT & current_estimate,
const int subset_num )
virtual

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

\[ {\partial L \over \partial \lambda_v} + P_v =
   \sum_b P_{bv} {y_b \over Y_b}
\]

(see the class general documentation). The sum will however be restricted to a subset.

Reimplemented in stir::PoissonLogLikelihoodWithLinearModelForMeanAndGatedProjDataWithMotion< TargetT >.

References actual_compute_subset_gradient_without_penalty(), compute_sub_gradient_without_penalty_plus_sensitivity(), and stir::error().

Referenced by compute_sub_gradient_without_penalty_plus_sensitivity().

◆ set_up()

template<typename TargetT>
Succeeded stir::PoissonLogLikelihoodWithLinearModelForMean< TargetT >::set_up ( shared_ptr< TargetT > const & target_sptr)
overridevirtual

set-up sensitivity etc if possible

If recompute_sensitivity is false, we will try to read it from either subset_sensitivity_filenames or sensitivity_filename, depending on the setting of get_use_subset_sensitivities().

If sensitivity_filename is equal to "1", all data are set to 1.

Deprecated
This special handling of the string "1" will be removed later.

Calls set_up_before_sensitivity().

Reimplemented from stir::GeneralisedObjectiveFunction< TargetT >.

References get_use_subset_sensitivities(), stir::GeneralisedObjectiveFunction< TargetT >::set_up(), and set_up().

Referenced by set_up().

◆ add_subset_sensitivity()

◆ get_use_subset_sensitivities()

template<typename TargetT>
bool stir::PoissonLogLikelihoodWithLinearModelForMean< TargetT >::get_use_subset_sensitivities ( ) const

find out if subset_sensitivities are used

If true, the sub_gradient and subset_sensitivity functions use the sensitivity for the given subset, otherwise, we use the total sensitivity divided by the number of subsets. The latter uses less memory, but is less stable for most (all?) algorithms.

References get_use_subset_sensitivities().

Referenced by stir::PoissonLogLikelihoodWithLinearModelForMeanAndListModeDataWithProjMatrixByBin< TargetT >::actual_compute_objective_function_without_penalty(), stir::PoissonLogLikelihoodWithLinearModelForMeanAndListModeDataWithProjMatrixByBin< TargetT >::actual_compute_subset_gradient_without_penalty(), add_subset_sensitivity(), compute_sensitivities(), get_use_subset_sensitivities(), and set_up().

◆ get_sensitivity_filename()

template<typename TargetT>
std::string stir::PoissonLogLikelihoodWithLinearModelForMean< TargetT >::get_sensitivity_filename ( ) const

get filename to read (or write) the total sensitivity

will be a zero-length string if not set

References get_sensitivity_filename().

Referenced by add_subset_sensitivity(), and get_sensitivity_filename().

◆ get_subsensitivity_filenames()

template<typename TargetT>
std::string stir::PoissonLogLikelihoodWithLinearModelForMean< TargetT >::get_subsensitivity_filenames ( ) const

get filename pattern to read (or write) the subset sensitivities

will be a zero-length string if not set. Could be e.g. "subsens_%d.hv" fmt::format is used with the pattern

References get_subsensitivity_filenames().

Referenced by add_subset_sensitivity(), and get_subsensitivity_filenames().

◆ get_subsensitivity_filename()

template<typename TargetT>
std::string stir::PoissonLogLikelihoodWithLinearModelForMean< TargetT >::get_subsensitivity_filename ( const int subset_num) const

◆ set_sensitivity_filename()

template<typename TargetT>
void stir::PoissonLogLikelihoodWithLinearModelForMean< TargetT >::set_sensitivity_filename ( const std::string & filename)

set filename to read (or write) the total sensitivity

set to a zero-length string to avoid reading/writing a file

References set_sensitivity_filename().

Referenced by add_subset_sensitivity(), and set_sensitivity_filename().

◆ set_subsensitivity_filenames()

template<typename TargetT>
void stir::PoissonLogLikelihoodWithLinearModelForMean< TargetT >::set_subsensitivity_filenames ( const std::string & filenames)

set filename pattern to read (or write) the subset sensitivities

set to a zero-length string to avoid reading/writing a file Could be e.g. "subsens_%d.hv" fmt::format is used with the pattern (which means you can use it like sprintf)

Calls error() if the pattern is invalid.

References stir::error(), get_subsensitivity_filename(), set_subsensitivity_filenames(), and stir::warning().

Referenced by add_subset_sensitivity(), post_processing(), stir::PoissonLogLikelihoodWithLinearModelForMeanAndListModeData< TargetT >::set_skip_lm_input_file(), and set_subsensitivity_filenames().

◆ fill_nonidentifiable_target_parameters()

template<typename TargetT>
void stir::PoissonLogLikelihoodWithLinearModelForMean< TargetT >::fill_nonidentifiable_target_parameters ( TargetT & target,
const float value ) const
overridevirtual

The implementation checks if the sensitivity of a voxel is zero. If so, it will the target voxel will be assigned the desired value.

Reimplemented from stir::GeneralisedObjectiveFunction< TargetT >.

References fill_nonidentifiable_target_parameters(), and get_sensitivity().

Referenced by add_subset_sensitivity(), and fill_nonidentifiable_target_parameters().

◆ set_up_before_sensitivity()

◆ compute_sensitivities()

template<typename TargetT>
void stir::PoissonLogLikelihoodWithLinearModelForMean< TargetT >::compute_sensitivities ( )
protected

compute subset and total sensitivity

This function fills in the sensitivity data by calling add_subset_sensitivity() for all subsets. It assumes that the subsensitivity for the 1st subset has been allocated already (and is the correct size).

References add_subset_sensitivity(), compute_sensitivities(), stir::error(), get_use_subset_sensitivities(), and stir::GeneralisedObjectiveFunction< TargetT >::subsets_are_approximately_balanced().

Referenced by compute_sensitivities(), and set_up_before_sensitivity().

◆ actual_compute_subset_gradient_without_penalty()

template<typename TargetT>
virtual void stir::PoissonLogLikelihoodWithLinearModelForMean< TargetT >::actual_compute_subset_gradient_without_penalty ( TargetT & gradient,
const TargetT & current_estimate,
const int subset_num,
const bool add_sensitivity )
protectedpure virtual

computes the subset gradient of the objective function without the penalty (optional: add subset sensitivity)

If add_sensitivity is true, this computes

\[ {\partial L \over \partial \lambda_v} + P_v =
  \sum_b P_{bv} {y_b \over Y_b}
\]

(see the class general documentation). The sum will however be restricted to a subset.

However, if add_sensitivity is false, this function will instead compute only the gradient

\[    {\partial L \over \partial \lambda_v} =
        \sum_b P_{bv} ({y_b \over Y_b} - 1)
\]

Implemented in stir::PoissonLogLikelihoodWithLinearKineticModelAndDynamicProjectionData< TargetT >, stir::PoissonLogLikelihoodWithLinearModelForMeanAndGatedProjDataWithMotion< TargetT >, stir::PoissonLogLikelihoodWithLinearModelForMeanAndListModeDataWithProjMatrixByBin< TargetT >, stir::PoissonLogLikelihoodWithLinearModelForMeanAndListModeDataWithProjMatrixByBin< stir::DiscretisedDensity< 3, float > >, and stir::PoissonLogLikelihoodWithLinearModelForMeanAndProjData< TargetT >.

References initialise_keymap(), post_processing(), and set_defaults().

Referenced by compute_sub_gradient_without_penalty(), and compute_sub_gradient_without_penalty_plus_sensitivity().

◆ set_defaults()

◆ initialise_keymap()

◆ post_processing()


The documentation for this class was generated from the following files: