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 >:

[legend]

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. More...

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. More...

Succeeded	set_up (shared_ptr< TargetT > const &target_sptr) override
	set-up sensitivity etc if possible More...

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 More...

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 More...

std::string	get_subsensitivity_filenames () const
	get filename pattern to read (or write) the subset sensitivities More...

void	fill_nonidentifiable_target_parameters (TargetT &target, const float value) const override

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 More...

void	set_subsensitivity_filenames (const std::string &)
	set filename pattern to read (or write) the subset sensitivities More...

Public Member Functions inherited from stir::GeneralisedObjectiveFunction< TargetT >
virtual TargetT *	construct_target_ptr () const =0
	Creates a suitable target as determined by the parameters. More...

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. More...

virtual void	compute_gradient (TargetT &gradient, const TargetT &current_estimate)
	Compute the gradient of the objective function at the current_estimate. More...

virtual void	compute_gradient_without_penalty (TargetT &gradient, const TargetT &current_estimate)
	Compute the gradient of the unregularised objective function at the current_estimate. More...

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. More...

virtual double	compute_objective_function_without_penalty (const TargetT &current_estimate)
	Compute the value of the unregularised objective function at the current_estimate. More...

double	compute_penalty (const TargetT &current_estimate, const int subset_num)
	Compute the value of the sub-penalty at the current_estimate. More...

double	compute_penalty (const TargetT &current_estimate)
	Compute the value of the penalty at the current_estimate. More...

double	compute_objective_function (const TargetT &current_estimate, const int subset_num)
	Compute the value of the sub-objective function at the current_estimate. More...

double	compute_objective_function (const TargetT &current_estimate)
	Compute the value of the objective function at the current_estimate. More...

double	compute_value (const TargetT &current_estimate)
	Alias for compute_objective_function(const TargetT&)

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< ExamInfo >	get_exam_info_uptr_for_target () const

virtual int	set_num_subsets (const int num_subsets)=0
	Attempts to change the number of subsets. More...

bool	subsets_are_approximately_balanced () const
	Checks of the current subset scheme is approximately balanced. More...

bool	subsets_are_approximately_balanced (std::string &warning_message) const
	Checks of the current subset scheme is approximately balanced and constructs a warning message. More...

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. More...

shared_ptr< GeneralisedPrior< TargetT > >	get_prior_sptr ()

void	set_prior_sptr (const shared_ptr< GeneralisedPrior< TargetT >> &)
	Change the prior. More...

virtual void	set_input_data (const shared_ptr< ExamData > &)=0
	set_input_data More...

virtual const ExamData &	get_input_data () const =0
	get input data More...

virtual void	set_additive_proj_data_sptr (const shared_ptr< ExamData > &)=0
	set_additive_proj_data_sptr More...

virtual void	set_normalisation_sptr (const shared_ptr< BinNormalisation > &)=0
	set_normalisation_sptr More...

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\'. More...

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\'. More...

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

Public Member Functions inherited from stir::RegisteredObjectBase
virtual std::string	get_registered_name () const =0
	Returns the name of the type of the object. More...

Public Member Functions inherited from stir::ParsingObject
	ParsingObject (const ParsingObject &)

ParsingObject &	operator= (const ParsingObject &)

void	ask_parameters ()

virtual std::string	parameter_info ()

bool	parse (std::istream &f)

bool	parse (const char *const filename)

Protected Member Functions
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 More...

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) More...

void	set_defaults () override
	Sets defaults for parsing. More...

void	initialise_keymap () override
	sets parsing keys More...

bool	post_processing () override
	This will be called at the end of the parsing. More...

Protected Member Functions inherited from stir::GeneralisedObjectiveFunction< TargetT >
virtual bool	actual_subsets_are_approximately_balanced (std::string &warning_message) const =0
	Implementation of function that checks subset balancing. More...

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. More...

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. More...

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. More...

Protected Member Functions inherited from stir::ParsingObject
virtual void	set_key_values ()
	This will be called before parsing or parameter_info is called. More...

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. More...

static GeneralisedObjectiveFunction< TargetT > *	ask_type_and_parameters ()
	ask the user for the type, and then calls read_registered_object(0, type) More...

static void	list_registered_names (std::ostream &stream)
	List all possible registered names to the stream. More...

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_less >	RegistryType
	The type of the registry.

Static Protected Member Functions inherited from stir::RegisteredObject< GeneralisedObjectiveFunction< TargetT > >
static RegistryType &	registry ()
	Static function returning the registry. More...

Protected Attributes inherited from stir::GeneralisedObjectiveFunction< TargetT >
int	num_subsets

bool	already_set_up

shared_ptr< GeneralisedPrior< TargetT > >	prior_sptr

Protected Attributes inherited from stir::ParsingObject
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.

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
; boost::format is used with the pattern (which means you can use it like sprintf)
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 >.

Reimplemented in stir::SumOfGeneralisedObjectiveFunctions< PoissonLogLikelihoodWithLinearModelForMeanAndProjData< TargetT >, TargetT, PoissonLogLikelihoodWithLinearModelForMean< TargetT > >.

References stir::error().

◆ 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 stir::error().

◆ 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 >.

Reimplemented in stir::SumOfGeneralisedObjectiveFunctions< PoissonLogLikelihoodWithLinearModelForMeanAndProjData< TargetT >, TargetT, PoissonLogLikelihoodWithLinearModelForMean< TargetT > >.

Referenced by stir::PoissonLogLikelihoodWithLinearModelForMeanAndListModeDataWithProjMatrixByBinTests::run_tests_for_objective_function(), and stir::PoissonLogLikelihoodWithLinearModelForMeanAndProjDataTests::test_approximate_Hessian_concavity().

◆ 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.

◆ 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 string if not set

◆ 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 string if not set. Could be e.g. "subsens_%d.hv" boost::format is used with the pattern (which means you can use it like sprintf)

◆ 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

◆ 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" boost::format is used with the pattern (which means you can use it like sprintf)

Calls error() if the pattern is invalid.

References stir::error().

◆ 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 >.

◆ 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 stir::error().

◆ 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::PoissonLogLikelihoodWithLinearModelForMeanAndProjData< TargetT >, stir::PoissonLogLikelihoodWithLinearKineticModelAndDynamicProjectionData< TargetT >, stir::PoissonLogLikelihoodWithLinearModelForMeanAndListModeDataWithProjMatrixByBin< TargetT >, and stir::PoissonLogLikelihoodWithLinearModelForMeanAndGatedProjDataWithMotion< TargetT >.

◆ set_defaults()

template<typename TargetT >

void stir::PoissonLogLikelihoodWithLinearModelForMean< TargetT >::set_defaults ( )

overrideprotectedvirtual

Sets defaults for parsing.

Resets sensitivity_filename, subset_sensitivity_filenames to empty, recompute_sensitivity to false, and use_subset_sensitivities to false.

Reimplemented from stir::GeneralisedObjectiveFunction< TargetT >.

Reimplemented in stir::PoissonLogLikelihoodWithLinearModelForMeanAndProjData< TargetT >, stir::PoissonLogLikelihoodWithLinearModelForMeanAndGatedProjDataWithMotion< TargetT >, stir::PoissonLogLikelihoodWithLinearModelForMeanAndListModeData< TargetT >, stir::PoissonLogLikelihoodWithLinearModelForMeanAndDynamicProjData< TargetT >, and stir::PoissonLogLikelihoodWithLinearModelForMeanAndListModeDataWithProjMatrixByBin< TargetT >.

◆ initialise_keymap()

template<typename TargetT >

void stir::PoissonLogLikelihoodWithLinearModelForMean< TargetT >::initialise_keymap ( )

overrideprotectedvirtual

sets parsing keys

Has to be called by initialise_keymap in the leaf-class

Reimplemented from stir::GeneralisedObjectiveFunction< TargetT >.

Reimplemented in stir::PoissonLogLikelihoodWithLinearModelForMeanAndProjData< TargetT >, stir::PoissonLogLikelihoodWithLinearModelForMeanAndGatedProjDataWithMotion< TargetT >, stir::PoissonLogLikelihoodWithLinearModelForMeanAndListModeData< TargetT >, stir::PoissonLogLikelihoodWithLinearModelForMeanAndDynamicProjData< TargetT >, and stir::PoissonLogLikelihoodWithLinearModelForMeanAndListModeDataWithProjMatrixByBin< TargetT >.

◆ post_processing()

template<typename TargetT >

bool stir::PoissonLogLikelihoodWithLinearModelForMean< TargetT >::post_processing ( )

overrideprotectedvirtual

This will be called at the end of the parsing.

Returns: false if everything OK, true if not

Reimplemented from stir::ParsingObject.

Reimplemented in stir::PoissonLogLikelihoodWithLinearModelForMeanAndProjData< TargetT >, stir::PoissonLogLikelihoodWithLinearModelForMeanAndGatedProjDataWithMotion< TargetT >, stir::PoissonLogLikelihoodWithLinearModelForMeanAndListModeData< TargetT >, stir::PoissonLogLikelihoodWithLinearModelForMeanAndDynamicProjData< TargetT >, and stir::PoissonLogLikelihoodWithLinearModelForMeanAndListModeDataWithProjMatrixByBin< TargetT >.

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

/home/sirfuser/devel/STIRdistrib/STIR/src/include/stir/KOSMAPOSL/KOSMAPOSLReconstruction.h
/home/sirfuser/devel/STIRdistrib/STIR/src/include/stir/recon_buildblock/PoissonLogLikelihoodWithLinearModelForMean.h
/home/sirfuser/devel/STIRdistrib/STIR/src/recon_buildblock/PoissonLogLikelihoodWithLinearModelForMean.cxx

Public Member Functions

Protected Member Functions

Additional Inherited Members

Detailed Description

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

Member Function Documentation

◆ compute_sub_gradient_without_penalty()

◆ compute_sub_gradient_without_penalty_plus_sensitivity()

◆ set_up()

◆ get_use_subset_sensitivities()

◆ get_sensitivity_filename()

◆ get_subsensitivity_filenames()

◆ set_sensitivity_filename()

◆ set_subsensitivity_filenames()

◆ fill_nonidentifiable_target_parameters()

◆ compute_sensitivities()

◆ actual_compute_subset_gradient_without_penalty()

◆ set_defaults()

◆ initialise_keymap()

◆ post_processing()

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