From dcf76a063ab230c76b0d7e9745a64f3bf0b61cd6 Mon Sep 17 00:00:00 2001 From: Kris Thielemans Date: Mon, 13 May 2024 14:10:24 +0100 Subject: [PATCH] LM obj-fun: add value() and tests - compute_objective_function() is now implemented as well - Tests are (almost) a duplicate of the ProjData version, but gradient-test is disabled as for this sample data it is too slow --- documentation/release_6.1.htm | 23 +- ...orMeanAndListModeDataWithProjMatrixByBin.h | 14 +- ...MeanAndListModeDataWithProjMatrixByBin.cxx | 68 ++++- src/recon_test/CMakeLists.txt | 4 + ...lForMeanAndListModeWithProjMatrixByBin.cxx | 283 ++++++++++++++++++ 5 files changed, 371 insertions(+), 21 deletions(-) create mode 100644 src/recon_test/test_PoissonLogLikelihoodWithLinearModelForMeanAndListModeWithProjMatrixByBin.cxx diff --git a/documentation/release_6.1.htm b/documentation/release_6.1.htm index 905568554d..f0512e672d 100644 --- a/documentation/release_6.1.htm +++ b/documentation/release_6.1.htm @@ -57,9 +57,15 @@

New functionality

Read TOF bin order from interfile header (see PR #1389)
  • - GeneralisedObjectiveFunction has 2 new papers to add compute full gradient + PoissonLogLikelihoodWithLinearModelForMeanAndListModeDataWithProjMatrixByBin can now compute the value + as well as accumulate_Hessian_times_input. +
    + PR #1418 +
    • +
  • + GeneralisedObjectiveFunction has 2 new members to compute the full gradient (compute_gradient and compute_gradient_without_penalty). Previously, - only subset gradient were available. + only subset gradients were available.
    PR #1418
    • @@ -68,6 +74,12 @@

    New functionality

    Changed functionality

    C++ tests

    + diff --git a/src/include/stir/recon_buildblock/PoissonLogLikelihoodWithLinearModelForMeanAndListModeDataWithProjMatrixByBin.h b/src/include/stir/recon_buildblock/PoissonLogLikelihoodWithLinearModelForMeanAndListModeDataWithProjMatrixByBin.h index 1bd095f150..be1950f229 100644 --- a/src/include/stir/recon_buildblock/PoissonLogLikelihoodWithLinearModelForMeanAndListModeDataWithProjMatrixByBin.h +++ b/src/include/stir/recon_buildblock/PoissonLogLikelihoodWithLinearModelForMeanAndListModeDataWithProjMatrixByBin.h @@ -74,6 +74,13 @@ class PoissonLogLikelihoodWithLinearModelForMeanAndListModeDataWithProjMatrixByB PoissonLogLikelihoodWithLinearModelForMeanAndListModeDataWithProjMatrixByBin(); + //! Computes the value of the objective function at the \a current_estimate. + /*! + \warning If add_sensitivity = false and use_subset_sensitivities = false will return an error + because the value will not be correct. Try use_subset_sensitivities = true. + */ + double actual_compute_objective_function_without_penalty(const TargetT& current_estimate, const int subset_num) override; + //! Computes the gradient of the objective function at the \a current_estimate overwriting \a gradient. /*! \warning If add_sensitivity = false and use_subset_sensitivities = false will return an error @@ -107,13 +114,6 @@ class PoissonLogLikelihoodWithLinearModelForMeanAndListModeDataWithProjMatrixByB #endif protected: - /*! \todo this function is not implemented yet and currently calls error() */ - double actual_compute_objective_function_without_penalty(const TargetT& current_estimate, const int subset_num) override - { // TODO - error("compute_objective_function_without_penalty Not implemented yet"); - return 0; - } - Succeeded set_up_before_sensitivity(shared_ptr const& target_sptr) override; void add_subset_sensitivity(TargetT& sensitivity, const int subset_num) const override; diff --git a/src/recon_buildblock/PoissonLogLikelihoodWithLinearModelForMeanAndListModeDataWithProjMatrixByBin.cxx b/src/recon_buildblock/PoissonLogLikelihoodWithLinearModelForMeanAndListModeDataWithProjMatrixByBin.cxx index 56c5ea2f9b..a10c68e1ca 100644 --- a/src/recon_buildblock/PoissonLogLikelihoodWithLinearModelForMeanAndListModeDataWithProjMatrixByBin.cxx +++ b/src/recon_buildblock/PoissonLogLikelihoodWithLinearModelForMeanAndListModeDataWithProjMatrixByBin.cxx @@ -751,13 +751,14 @@ constexpr float max_quotient = 10000.F; \sum_e A_e^t (y_e/(A_e lambda+ c)) */ +template inline void -LM_gradient(DiscretisedDensity<3, float>& output_image, - const ProjMatrixElemsForOneBin& row, - const float add_term, - const Bin& measured_bin, - const DiscretisedDensity<3, float>& input_image, - double* value_ptr) +LM_gradient_and_value(DiscretisedDensity<3, float>& output_image, + const ProjMatrixElemsForOneBin& row, + const float add_term, + const Bin& measured_bin, + const DiscretisedDensity<3, float>& input_image, + double* value_ptr) { Bin fwd_bin = measured_bin; fwd_bin.set_bin_value(0.0f); @@ -767,18 +768,22 @@ LM_gradient(DiscretisedDensity<3, float>& output_image, if (measured_bin.get_bin_value() > max_quotient * fwd) { // cancel singularity - if (value_ptr) + if (do_value) { + assert(value_ptr); const auto num = measured_bin.get_bin_value(); *value_ptr -= num * log(double(num / max_quotient)); return; } } - const auto measured_div_fwd = measured_bin.get_bin_value() / fwd; + if (do_gradient) + { + const auto measured_div_fwd = measured_bin.get_bin_value() / fwd; - fwd_bin.set_bin_value(measured_div_fwd); - row.back_project(output_image, fwd_bin); - if (value_ptr) + fwd_bin.set_bin_value(measured_div_fwd); + row.back_project(output_image, fwd_bin); + } + if (do_value) *value_ptr -= measured_bin.get_bin_value() * log(double(fwd)); } @@ -835,7 +840,7 @@ LM_gradient_distributable_computation(const shared_ptr PM_sptr, has_add, accumulate, value_ptr, - LM_gradient); + LM_gradient_and_value); } void @@ -865,6 +870,45 @@ LM_Hessian_distributable_computation(const shared_ptr PM_sptr, H_func); } +template +double +PoissonLogLikelihoodWithLinearModelForMeanAndListModeDataWithProjMatrixByBin< + TargetT>::actual_compute_objective_function_without_penalty(const TargetT& current_estimate, const int subset_num) +{ + assert(subset_num >= 0); + assert(subset_num < this->num_subsets); + if (!this->get_use_subset_sensitivities() && this->num_subsets > 1) + error("PoissonLogLikelihoodWithLinearModelForMeanAndListModeDataWithProjMatrixByBin::" + "actual_compute_subset_gradient_without_penalty(): cannot subtract subset sensitivity because " + "use_subset_sensitivities is false. This will result in an error in the gradient computation."); + + double accum = 0.; + unsigned int icache = 0; + while (true) + { + bool stop = this->load_listmode_batch(icache); + LM_distributable_computation(this->PM_sptr, + this->proj_data_info_sptr, + nullptr, + ¤t_estimate, + record_cache, + subset_num, + this->num_subsets, + this->has_add, + /* accumulate */ true, + &accum, + LM_gradient_and_value); + ++icache; + if (stop) + break; + } + std::inner_product(current_estimate.begin_all_const(), + current_estimate.end_all_const(), + this->get_subset_sensitivity(subset_num).begin_all_const(), + accum); + return accum; +} + template void PoissonLogLikelihoodWithLinearModelForMeanAndListModeDataWithProjMatrixByBin< diff --git a/src/recon_test/CMakeLists.txt b/src/recon_test/CMakeLists.txt index f7fe116313..6e9e329f13 100644 --- a/src/recon_test/CMakeLists.txt +++ b/src/recon_test/CMakeLists.txt @@ -36,6 +36,7 @@ set(${dir_INVOLVED_TEST_EXE_SOURCES} recontest.cxx test_data_processor_projectors.cxx test_OSMAPOSL.cxx + test_PoissonLogLikelihoodWithLinearModelForMeanAndListModeWithProjMatrixByBin.cxx ) if (DOWNLOAD_ZENODO_TEST_DATA) @@ -128,6 +129,9 @@ include(stir_test_exe_targets) # a test that uses MPI create_stir_mpi_test(test_PoissonLogLikelihoodWithLinearModelForMeanAndProjData.cxx "${STIR_LIBRARIES}" $) +ADD_TEST(test_PoissonLogLikelihoodWithLinearModelForMeanAndListModeWithProjMatrixByBin + test_PoissonLogLikelihoodWithLinearModelForMeanAndListModeWithProjMatrixByBin "${CMAKE_SOURCE_DIR}/recon_test_pack/PET_ACQ_small.l.hdr.STIR") + # fwdtest and bcktest could be useful on their own, so we'll add them to the installation targets if (BUILD_TESTING) install(TARGETS fwdtest bcktest DESTINATION bin) diff --git a/src/recon_test/test_PoissonLogLikelihoodWithLinearModelForMeanAndListModeWithProjMatrixByBin.cxx b/src/recon_test/test_PoissonLogLikelihoodWithLinearModelForMeanAndListModeWithProjMatrixByBin.cxx new file mode 100644 index 0000000000..a8c6ba669c --- /dev/null +++ b/src/recon_test/test_PoissonLogLikelihoodWithLinearModelForMeanAndListModeWithProjMatrixByBin.cxx @@ -0,0 +1,283 @@ +/* + Copyright (C) 2011, Hammersmith Imanet Ltd + Copyright (C) 2013, 2021, 2024 University College London + This file is part of STIR. + + SPDX-License-Identifier: Apache-2.0 + + See STIR/LICENSE.txt for details +*/ +/*! + + \file + \ingroup recon_test + + \brief Test program for stir::PoissonLogLikelihoodWithLinearModelForMeanAndListModeDataWithProjMatrixByBin + + \par Usage + + 
    +  test_PoissonLogLikelihoodWithLinearModelForMeanAndListModeDataWithProjMatrixByBin lm_data_filename [ density_filename ]
+
    + where the last arguments are optional. See the class documentation for more info. + + \author Kris Thielemans + \author Robert Twyman Skelly +*/ + +#include "stir/recon_buildblock/test/ObjectiveFunctionTests.h" +#include "stir/VoxelsOnCartesianGrid.h" +#include "stir/listmode/CListModeData.h" +#include "stir/ExamInfo.h" +#include "stir/ProjDataInfo.h" +#include "stir/ProjDataInMemory.h" +#include "stir/SegmentByView.h" +#include "stir/Scanner.h" +#include "stir/recon_buildblock/PoissonLogLikelihoodWithLinearModelForMeanAndListModeDataWithProjMatrixByBin.h" +#include "stir/recon_buildblock/ProjMatrixByBinUsingRayTracing.h" +#include "stir/recon_buildblock/ProjectorByBinPairUsingProjMatrixByBin.h" +#include "stir/recon_buildblock/BinNormalisationFromProjData.h" +#include "stir/recon_buildblock/TrivialBinNormalisation.h" +//#include "stir/OSMAPOSL/OSMAPOSLReconstruction.h" +#include "stir/recon_buildblock/distributable_main.h" +#include "stir/IO/read_from_file.h" +#include "stir/IO/write_to_file.h" +#include "stir/info.h" +#include "stir/Succeeded.h" +#include "stir/num_threads.h" +#include +#include +#include +#include +#include +#include + +START_NAMESPACE_STIR + +/*! + \ingroup test + \brief Test class for PoissonLogLikelihoodWithLinearModelForMeanAndListModeDataWithProjMatrixByBin + + This is a somewhat preliminary implementation of a test that compares the result + of GeneralisedObjectiveFunction::compute_gradient + with a numerical gradient computed by using the + GeneralisedObjectiveFunction::compute_objective_function() function. + + The trouble with this is that compute the gradient voxel by voxel is obviously + terribly slow. A solution (for the test) would be to compute it only in + a subset of voxels or so. We'll leave this for later. + + Note that the test only works if the objective function is well-defined. For example, + if certain projections are non-zero, while the model estimates them to be zero, the + Poisson objective function is in theory infinite. + PoissonLogLikelihoodWithLinearModelForMeanAndListModeDataWithProjMatrixByBin uses some thresholds to try to + avoid overflow, but if there are too many of these bins, the total objective + function will become infinite. The numerical gradient then becomes ill-defined + (even in voxels that do not contribute to these bins). + +*/ +class PoissonLogLikelihoodWithLinearModelForMeanAndListModeDataWithProjMatrixByBinTests + : public ObjectiveFunctionTests< + PoissonLogLikelihoodWithLinearModelForMeanAndListModeDataWithProjMatrixByBin>, + DiscretisedDensity<3, float>> +{ +public: + //! Constructor that can take some input data to run the test with + /*! This makes it possible to run the test with your own data. However, beware that + it is very easy to set up a very long computation. See also the note about + non-zero measured bins. + + \todo it would be better to parse an objective function. That would allow us to set + all parameters from the command line. + */ + PoissonLogLikelihoodWithLinearModelForMeanAndListModeDataWithProjMatrixByBinTests(char const* const lm_data_filename, + char const* const density_filename = 0); + void construct_input_data(shared_ptr& density_sptr); + + void run_tests() override; + +protected: + char const* lm_data_filename; + char const* density_filename; + shared_ptr lm_data_sptr; + shared_ptr mult_proj_data_sptr; + shared_ptr add_proj_data_sptr; + shared_ptr> objective_function_sptr; + + //! run the test + void run_tests_for_objective_function(objective_function_type& objective_function, target_type& target); +}; + +PoissonLogLikelihoodWithLinearModelForMeanAndListModeDataWithProjMatrixByBinTests:: + PoissonLogLikelihoodWithLinearModelForMeanAndListModeDataWithProjMatrixByBinTests(char const* lm_data_filename, + char const* const density_filename) + : lm_data_filename(lm_data_filename), + density_filename(density_filename) +{} + +void +PoissonLogLikelihoodWithLinearModelForMeanAndListModeDataWithProjMatrixByBinTests::run_tests_for_objective_function( + objective_function_type& objective_function, target_type& target) +{ + // TODO enable this, but far too slow ATM + // std::cerr << "----- testing Gradient\n"; + // test_gradient("PoissonLLProjData", objective_function, target, 0.01F); + + std::cerr << "----- testing concavity via Hessian-vector product (accumulate_Hessian_times_input)\n"; + test_Hessian_concavity("PoissonLLProjData", objective_function, target); + + std::cerr << "----- testing Hessian-vector product (accumulate_Hessian_times_input)\n"; + test_Hessian("PoissonLLProjData", objective_function, target, 0.5F); +} + +void +PoissonLogLikelihoodWithLinearModelForMeanAndListModeDataWithProjMatrixByBinTests::construct_input_data( + shared_ptr& density_sptr) +{ + lm_data_sptr = read_from_file(this->lm_data_filename); + + if (this->density_filename == 0) + { + // construct a small image + + CartesianCoordinate3D origin(0, 0, 0); + const float zoom = 0.03F; + + density_sptr.reset(new VoxelsOnCartesianGrid( + lm_data_sptr->get_exam_info_sptr(), *lm_data_sptr->get_proj_data_info_sptr(), zoom, origin)); + write_to_file("target.hv", *density_sptr); + // fill with random numbers between 0 and 1 + typedef boost::mt19937 base_generator_type; + // initialize by reproducible seed + static base_generator_type generator(boost::uint32_t(42)); + static boost::uniform_01 random01(generator); + for (target_type::full_iterator iter = density_sptr->begin_all(); iter != density_sptr->end_all(); ++iter) + *iter = static_cast(random01()); + } + else + { + shared_ptr aptr(read_from_file(this->density_filename)); + density_sptr = aptr; + } + + // make odd to avoid difficulties with outer-bin that isn't filled-in when using symmetries + { + BasicCoordinate<3, int> min_ind, max_ind; + if (density_sptr->get_regular_range(min_ind, max_ind)) + { + for (int d = 2; d <= 3; ++d) + { + min_ind[d] = std::min(min_ind[d], -max_ind[d]); + max_ind[d] = std::max(-min_ind[d], max_ind[d]); + } + density_sptr->grow(IndexRange<3>(min_ind, max_ind)); + } + } + + auto proj_data_info_sptr = lm_data_sptr->get_proj_data_info_sptr()->create_shared_clone(); + // multiplicative term + shared_ptr bin_norm_sptr(new TrivialBinNormalisation()); + { + + mult_proj_data_sptr.reset(new ProjDataInMemory(lm_data_sptr->get_exam_info_sptr(), proj_data_info_sptr)); + for (int seg_num = proj_data_info_sptr->get_min_segment_num(); seg_num <= proj_data_info_sptr->get_max_segment_num(); + ++seg_num) + { + for (int timing_pos_num = proj_data_info_sptr->get_min_tof_pos_num(); + timing_pos_num <= proj_data_info_sptr->get_max_tof_pos_num(); + ++timing_pos_num) + { + SegmentByView segment = proj_data_info_sptr->get_empty_segment_by_view(seg_num, false, timing_pos_num); + // fill in some crazy values + float value = 0; + for (SegmentByView::full_iterator iter = segment.begin_all(); iter != segment.end_all(); ++iter) + { + value = float(fabs(seg_num * value - .2)); // needs to be positive for Poisson + *iter = value; + } + segment /= 0.5F * segment.find_max(); // normalise to 2 (to avoid large numbers) + mult_proj_data_sptr->set_segment(segment); + } + } + bin_norm_sptr.reset(new BinNormalisationFromProjData(mult_proj_data_sptr)); + } + + // additive term + add_proj_data_sptr.reset(new ProjDataInMemory(lm_data_sptr->get_exam_info_sptr(), proj_data_info_sptr)); + { + for (int seg_num = proj_data_info_sptr->get_min_segment_num(); seg_num <= proj_data_info_sptr->get_max_segment_num(); + ++seg_num) + { + for (int timing_pos_num = proj_data_info_sptr->get_min_tof_pos_num(); + timing_pos_num <= proj_data_info_sptr->get_max_tof_pos_num(); + ++timing_pos_num) + { + SegmentByView segment = proj_data_info_sptr->get_empty_segment_by_view(seg_num, false, timing_pos_num); + // fill in some crazy values + float value = 0; + for (SegmentByView::full_iterator iter = segment.begin_all(); iter != segment.end_all(); ++iter) + { + value = float(fabs(seg_num * value - .3)); // needs to be positive for Poisson + *iter = value; + } + segment /= 0.333F * segment.find_max(); // normalise to 3 (to avoid large numbers) + add_proj_data_sptr->set_segment(segment); + } + } + } + + objective_function_sptr.reset(new PoissonLogLikelihoodWithLinearModelForMeanAndListModeDataWithProjMatrixByBin); + PoissonLogLikelihoodWithLinearModelForMeanAndListModeDataWithProjMatrixByBin& objective_function + = reinterpret_cast&>( + *objective_function_sptr); + objective_function.set_input_data(lm_data_sptr); + objective_function.set_use_subset_sensitivities(true); + objective_function.set_max_segment_num_to_process(1); + shared_ptr proj_matrix_sptr(new ProjMatrixByBinUsingRayTracing()); + objective_function.set_proj_matrix(proj_matrix_sptr); + objective_function.set_normalisation_sptr(bin_norm_sptr); + objective_function.set_additive_proj_data_sptr(add_proj_data_sptr); + objective_function.set_num_subsets(2); + if (!check(objective_function.set_up(density_sptr) == Succeeded::yes, "set-up of objective function")) + return; +} + +void +PoissonLogLikelihoodWithLinearModelForMeanAndListModeDataWithProjMatrixByBinTests::run_tests() +{ + std::cerr << "Tests for PoissonLogLikelihoodWithLinearModelForMeanAndListModeDataWithProjMatrixByBin\n"; + const int verbosity_default = Verbosity::get(); + Verbosity::set(2); + +#if 1 + shared_ptr density_sptr; + construct_input_data(density_sptr); + this->run_tests_for_objective_function(*this->objective_function_sptr, *density_sptr); +#else + // alternative that gets the objective function from an OSMAPOSL .par file + // currently disabled + OSMAPOSLReconstruction recon(proj_data_filename); // actually .par + shared_ptr> objective_function_sptr = recon.get_objective_function_sptr(); + if (!check(objective_function_sptr->set_up(recon.get_initial_data_ptr()) == Succeeded::yes, "set-up of objective function")) + return; + this->run_tests_for_objective_function(*objective_function_sptr, *recon.get_initial_data_ptr()); +#endif + Verbosity::set(verbosity_default); +} + +END_NAMESPACE_STIR + +USING_NAMESPACE_STIR + +int +main(int argc, char** argv) +{ + if (argc <= 1) + error("Need to specify a list-mode filename"); + + set_default_num_threads(); + + PoissonLogLikelihoodWithLinearModelForMeanAndListModeDataWithProjMatrixByBinTests tests(argv[1], argc > 2 ? argv[2] : 0); + tests.run_tests(); + return tests.main_return_value(); +}