From 267c16123570b1ef995e634c8554f42b966966f1 Mon Sep 17 00:00:00 2001 From: Kris Thielemans Date: Tue, 14 May 2024 08:36:09 +0100 Subject: [PATCH] add "increment" test for gradient (enabled for LM data) also minor clean-up --- .../test/ObjectiveFunctionTests.h | 105 ++++++++++++------ src/recon_test/CMakeLists.txt | 1 + ...lForMeanAndListModeWithProjMatrixByBin.cxx | 9 +- ...ihoodWithLinearModelForMeanAndProjData.cxx | 6 +- 4 files changed, 79 insertions(+), 42 deletions(-) diff --git a/src/include/stir/recon_buildblock/test/ObjectiveFunctionTests.h b/src/include/stir/recon_buildblock/test/ObjectiveFunctionTests.h index 58a53731c..fce90ef36 100644 --- a/src/include/stir/recon_buildblock/test/ObjectiveFunctionTests.h +++ b/src/include/stir/recon_buildblock/test/ObjectiveFunctionTests.h @@ -15,7 +15,7 @@ \brief Test skeleton for stir::GeneralisedObjectiveFunction and stir::GeneralisedPrior \author Kris Thielemans - \author Reobert Twyman Skelly + \author Robert Twyman Skelly */ #include "stir/RunTests.h" @@ -32,10 +32,6 @@ START_NAMESPACE_STIR This contains some numerical tests to check gradient and Hessian calculations. - Note that the gradient is computed numerically voxel by voxel which 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. @@ -53,16 +49,22 @@ class ObjectiveFunctionTests : public RunTests typedef TargetT target_type; //! Test the gradient of the objective function by comparing to the numerical gradient via perturbation - /*! Note: \a target is non-\c const, as the code will add/subtract eps, but the actual values + /*! + If \a full_gradient=true, all elements in the gradient are tested (using single-element increments). This is slow. + Otherwise, the test checks that \f$ G dx \approx F(x+dx) - F(x) \f$. dx is computed via construct_increment(). + + Note: \a target is non-\c const, as the code will add/subtract eps, but the actual values are not modified after the test exits. */ - virtual Succeeded - test_gradient(const std::string& test_name, ObjectiveFunctionT& objective_function, TargetT& target, const float eps); + virtual Succeeded test_gradient(const std::string& test_name, + ObjectiveFunctionT& objective_function, + TargetT& target, + const float eps, + const bool full_gradient = true); //! Test the accumulate_Hessian_times_input of the objective function by comparing to the numerical result via perturbation /*! - This test checks that \f$ H dx \approx G(x+dx) - G(x) \f$. - \f$dx\f$ is currently computed by as eps*(target / target.find_max() + 0.5). + This test checks that \f$ H dx \approx G(x+dx) - G(x) \f$. dx is computed via construct_increment(). */ virtual Succeeded test_Hessian(const std::string& test_name, ObjectiveFunctionT& objective_function, const TargetT& target, const float eps); @@ -71,6 +73,14 @@ class ObjectiveFunctionTests : public RunTests ObjectiveFunctionT& objective_function, const TargetT& target, const float mult_factor = 1.F); + +protected: + //! Construct small increment for target + /*! + Result is eps*(target / target.find_max() + 0.5), i.e. it is always + positive (if target is non-negative). + */ + virtual shared_ptr construct_increment(const TargetT& target, const float eps) const; }; template @@ -78,7 +88,8 @@ Succeeded ObjectiveFunctionTests::test_gradient(const std::string& test_name, ObjectiveFunctionT& objective_function, TargetT& target, - const float eps) + const float eps, + const bool full_gradient) { shared_ptr gradient_sptr(target.get_empty_copy()); shared_ptr gradient_2_sptr(target.get_empty_copy()); @@ -87,29 +98,47 @@ ObjectiveFunctionTests::test_gradient(const std::st this->set_tolerance(std::max(fabs(double(gradient_sptr->find_min())), double(gradient_sptr->find_max())) / 1000); info("Computing objective function at target"); const double value_at_target = objective_function.compute_value(target); - auto target_iter = target.begin_all(); - auto gradient_iter = gradient_sptr->begin_all(); - auto gradient_2_iter = gradient_2_sptr->begin_all(); bool testOK = true; - info("Computing gradient of objective function by numerical differences (this will take a while)"); - while (target_iter != target.end_all()) + if (full_gradient) { - *target_iter += eps; - const double value_at_inc = objective_function.compute_value(target); - *target_iter -= eps; - const float gradient_at_iter = static_cast((value_at_inc - value_at_target) / eps); - *gradient_2_iter++ = gradient_at_iter; - testOK = testOK && this->check_if_equal(gradient_at_iter, *gradient_iter, "gradient"); - ++target_iter; - ++gradient_iter; + info("Computing gradient of objective function by numerical differences (this will take a while)"); + auto target_iter = target.begin_all(); + auto gradient_iter = gradient_sptr->begin_all(); + auto gradient_2_iter = gradient_2_sptr->begin_all(); + while (target_iter != target.end_all()) + { + *target_iter += eps; + const double value_at_inc = objective_function.compute_value(target); + *target_iter -= eps; + const float gradient_at_iter = static_cast((value_at_inc - value_at_target) / eps); + *gradient_2_iter++ = gradient_at_iter; + testOK = testOK && this->check_if_equal(gradient_at_iter, *gradient_iter, "gradient"); + ++target_iter; + ++gradient_iter; + } + } + else + { + /* test f(x+dx) - f(x) ~ dx^t G(x) */ + shared_ptr increment_sptr = this->construct_increment(target, eps); + shared_ptr target_plus_inc_sptr(target.clone()); + *target_plus_inc_sptr += *increment_sptr; + const double value_at_inc = objective_function.compute_value(*target_plus_inc_sptr); + const double my_sum = std::inner_product( + gradient_sptr->begin_all_const(), gradient_sptr->end_all_const(), increment_sptr->begin_all_const(), 0.); + + testOK = testOK && this->check_if_equal(value_at_inc - value_at_target, my_sum, "gradient"); } + if (!testOK) { std::cerr << "Numerical gradient test failed with for " + test_name + "\n"; - std::cerr << "Writing diagnostic files " << test_name << "_target.hv, *gradient.hv, *numerical_gradient.hv\n"; + std::cerr << "Writing diagnostic files " << test_name + << "_target.hv, *gradient.hv (and *numerical_gradient.hv if full gradient test is used)\n"; write_to_file(test_name + "_target.hv", target); write_to_file(test_name + "_gradient.hv", *gradient_sptr); - write_to_file(test_name + "_numerical_gradient.hv", *gradient_2_sptr); + if (full_gradient) + write_to_file(test_name + "_numerical_gradient.hv", *gradient_2_sptr); return Succeeded::no; } else @@ -118,6 +147,16 @@ ObjectiveFunctionTests::test_gradient(const std::st } } +template +shared_ptr +ObjectiveFunctionTests::construct_increment(const TargetT& target, const float eps) const +{ + shared_ptr increment_sptr(target.clone()); + *increment_sptr *= eps / increment_sptr->find_max(); + *increment_sptr += eps / 2; + return increment_sptr; +} + template Succeeded ObjectiveFunctionTests::test_Hessian(const std::string& test_name, @@ -125,13 +164,13 @@ ObjectiveFunctionTests::test_Hessian(const std::str const TargetT& target, const float eps) { + info("Comparing Hessian*dx with difference of gradients"); + /* test G(x+dx) = G(x) + H dx + small stuff */ shared_ptr gradient_sptr(target.get_empty_copy()); shared_ptr gradient_2_sptr(target.get_empty_copy()); shared_ptr output(target.get_empty_copy()); - shared_ptr increment_sptr(target.clone()); - *increment_sptr *= eps / increment_sptr->find_max(); - *increment_sptr += eps / 2; + shared_ptr increment_sptr = this->construct_increment(target, eps); shared_ptr target_plus_inc_sptr(target.clone()); *target_plus_inc_sptr += *increment_sptr; @@ -145,7 +184,6 @@ ObjectiveFunctionTests::test_Hessian(const std::str auto gradient_iter = gradient_sptr->begin_all_const(); auto gradient_2_iter = gradient_2_sptr->begin_all_const(); bool testOK = true; - info("Computing gradient of objective function by numerical differences (this will take a while)"); while (output_iter != output->end_all()) { testOK = testOK && this->check_if_equal(*gradient_2_iter - *gradient_iter, *output_iter, "Hessian*increment"); @@ -187,10 +225,9 @@ ObjectiveFunctionTests::test_Hessian_concavity(cons /// Compute dot(x,(H x)) const float my_sum = std::inner_product(target.begin_all(), target.end_all(), output->begin_all(), 0.F) * mult_factor; - // test for a CONCAVE function + // test for a CONCAVE function (after multiplying with mult_factor of course) if (this->check_if_less(my_sum, 0)) { - // info("PASS: Computation of x^T H x = " + std::to_string(my_sum) + " < 0" (Hessian) and is therefore concave); return Succeeded::yes; } else @@ -199,8 +236,8 @@ ObjectiveFunctionTests::test_Hessian_concavity(cons std::cerr << "FAIL: " + test_name + ": Computation of x^T H x = " + std::to_string(my_sum) << " > 0 (Hessian) and is therefore NOT concave" << "\n >target image max=" << target.find_max() << "\n >target image min=" << target.find_min() - << "\n >output image max=" << output->find_max() << "\n >output image min=" << output->find_min(); - std::cerr << "Writing diagnostic files to " << test_name + "_concavity_out.hv etc\n."; + << "\n >output image max=" << output->find_max() << "\n >output image min=" << output->find_min() << '\n'; + std::cerr << "Writing diagnostic files to " << test_name + "_concavity_out.hv, *target.hv\n"; write_to_file(test_name + "_concavity_out.hv", *output); write_to_file(test_name + "_target.hv", target); return Succeeded::no; diff --git a/src/recon_test/CMakeLists.txt b/src/recon_test/CMakeLists.txt index 6e9e329f1..8a5eaf169 100644 --- a/src/recon_test/CMakeLists.txt +++ b/src/recon_test/CMakeLists.txt @@ -129,6 +129,7 @@ include(stir_test_exe_targets) # a test that uses MPI create_stir_mpi_test(test_PoissonLogLikelihoodWithLinearModelForMeanAndProjData.cxx "${STIR_LIBRARIES}" $) +# pass list-mode file as argument ADD_TEST(test_PoissonLogLikelihoodWithLinearModelForMeanAndListModeWithProjMatrixByBin test_PoissonLogLikelihoodWithLinearModelForMeanAndListModeWithProjMatrixByBin "${CMAKE_SOURCE_DIR}/recon_test_pack/PET_ACQ_small.l.hdr.STIR") diff --git a/src/recon_test/test_PoissonLogLikelihoodWithLinearModelForMeanAndListModeWithProjMatrixByBin.cxx b/src/recon_test/test_PoissonLogLikelihoodWithLinearModelForMeanAndListModeWithProjMatrixByBin.cxx index a8c6ba669..46b2bce25 100644 --- a/src/recon_test/test_PoissonLogLikelihoodWithLinearModelForMeanAndListModeWithProjMatrixByBin.cxx +++ b/src/recon_test/test_PoissonLogLikelihoodWithLinearModelForMeanAndListModeWithProjMatrixByBin.cxx @@ -119,15 +119,14 @@ 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 Gradient\n"; + test_gradient("PoissonLLListModeData", objective_function, target, 0.01F, /* full_gradient = */ false); std::cerr << "----- testing concavity via Hessian-vector product (accumulate_Hessian_times_input)\n"; - test_Hessian_concavity("PoissonLLProjData", objective_function, target); + test_Hessian_concavity("PoissonLLListModeData", objective_function, target); std::cerr << "----- testing Hessian-vector product (accumulate_Hessian_times_input)\n"; - test_Hessian("PoissonLLProjData", objective_function, target, 0.5F); + test_Hessian("PoissonLLListModeData", objective_function, target, 0.5F); } void diff --git a/src/recon_test/test_PoissonLogLikelihoodWithLinearModelForMeanAndProjData.cxx b/src/recon_test/test_PoissonLogLikelihoodWithLinearModelForMeanAndProjData.cxx index 8fb8ee4a7..73d47fece 100644 --- a/src/recon_test/test_PoissonLogLikelihoodWithLinearModelForMeanAndProjData.cxx +++ b/src/recon_test/test_PoissonLogLikelihoodWithLinearModelForMeanAndProjData.cxx @@ -124,16 +124,16 @@ PoissonLogLikelihoodWithLinearModelForMeanAndProjDataTests::run_tests_for_object objective_function_type& objective_function, target_type& target) { std::cerr << "----- testing Gradient\n"; - test_gradient("PoissonLLListModeData", objective_function, target, 0.01F); + test_gradient("PoissonLLProjData", objective_function, target, 0.01F); std::cerr << "----- testing concavity via Hessian-vector product (accumulate_Hessian_times_input)\n"; - test_Hessian_concavity("PoissonLLListModeData", objective_function, target); + test_Hessian_concavity("PoissonLLProjData", objective_function, target); std::cerr << "----- testing approximate-Hessian-vector product (accumulate_Hessian_times_input)\n"; test_approximate_Hessian_concavity(objective_function, target); std::cerr << "----- testing Hessian-vector product (accumulate_Hessian_times_input)\n"; - test_Hessian("PoissonLLListModeData", objective_function, target, 0.5F); + test_Hessian("PoissonLLProjData", objective_function, target, 0.5F); if (!this->is_everything_ok()) {