From 7031363d390befd8f024a593c9365f62760a2de6 Mon Sep 17 00:00:00 2001 From: Jürgen Knödlseder Date: Fri, 25 Jun 2021 14:01:54 +0200 Subject: [PATCH] Add restart of bisection in case of bad convergence in ctulimit --- ChangeLog | 2 +- NEWS | 9 +++++---- doc/source/admin/release_history/2.0.rst | 3 ++- src/ctulimit/ctulimit.cpp | 101 +++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++------------------------------ src/ctulimit/ctulimit.hpp | 4 +++- 5 files changed, 82 insertions(+), 37 deletions(-) diff --git a/ChangeLog b/ChangeLog index 790690a..bb7244f 100755 --- a/ChangeLog +++ b/ChangeLog @@ -3,7 +3,7 @@ * Version 2.0.0 released ======================== - Re-compute best log-likelihood in ctulimit. + Add restart of bisection in case of bad convergence in ctulimit Add COMPTEL analysis scripts (#3738) Add source counts constraints to cssens (#3668) Stamp all output FITS files (#3662) diff --git a/NEWS b/NEWS index 9ac67ec..e951c5e 100755 --- a/NEWS +++ b/NEWS @@ -155,10 +155,11 @@ model parameters (#3562). ctulimit - Compute upper limit ------------------------------ -The best log-likelihood is recomputed for a fixed best parameter value to -make sure that the value is compliant with the values obtained by the bisection -algorithm. The bisection iterations are stopped with a warning if the parameter -interval becomes tiny. +If the bisection results in a result that is incompatible with the previously +determine best maximum log-likelihood value or if the bisection interval tends +towards zero then update the best log-likelihood value with the current result +and restart the bisection algorithm once more. This allows for dealing with +ill-conditioned fit problems. support - Support classes diff --git a/doc/source/admin/release_history/2.0.rst b/doc/source/admin/release_history/2.0.rst index c1636c5..7050be9 100644 --- a/doc/source/admin/release_history/2.0.rst +++ b/doc/source/admin/release_history/2.0.rst @@ -1,6 +1,6 @@ .. _2.0: -ctools 2.0 release (24 June 2021) +ctools 2.0 release (25 June 2021) ================================= Overview @@ -24,6 +24,7 @@ Bug fixes Improvements ------------ +* Add restart of bisection in case of bad convergence in :ref:`ctulimit` * [`3662 `_] - Stamp all output FITS files * [`3668 `_] - diff --git a/src/ctulimit/ctulimit.cpp b/src/ctulimit/ctulimit.cpp index 725a035..2945e6b 100755 --- a/src/ctulimit/ctulimit.cpp +++ b/src/ctulimit/ctulimit.cpp @@ -280,17 +280,18 @@ void ctulimit::run(void) // Write some parameters into logger log_value(NORMAL, "Model name", m_skymodel->name()); log_value(NORMAL, "Parameter name", m_model_par->name()); + log_value(NORMAL, "Parameter value", m_model_par->value()); log_value(NORMAL, "Confidence level", gammalib::str(m_confidence*100.0)+"%"); - log_value(NORMAL, "Best log-likelihood value", m_best_logL); log_value(NORMAL, "Log-likelihood difference", m_dlogL); - log_value(NORMAL, "Initial parameter range", + log_value(NORMAL, "Initial factor value range", "["+gammalib::str(parmin)+", "+gammalib::str(parmax)+"]"); // Compute upper limit ulimit_bisection(parmin, parmax); // Write final parameter into logger - log_value(NORMAL, "Final parameter", m_model_par->value()); + log_value(NORMAL, "Upper limit factor value", m_model_par->factor_value()); + log_value(NORMAL, "Upper limit param. value", m_model_par->value()); // Compute upper limit intensity and fluxes compute_ulimit(); @@ -370,6 +371,8 @@ void ctulimit::init_members(void) m_skymodel = NULL; m_model_par = NULL; m_best_logL = 0.0; + m_best_value = 0.0; + m_best_error = 0.0; m_flux_ulimit = 0.0; m_diff_ulimit = 0.0; m_eflux_ulimit = 0.0; @@ -403,6 +406,8 @@ void ctulimit::copy_members(const ctulimit& app) // Copy protected members m_dlogL = app.m_dlogL; m_best_logL = app.m_best_logL; + m_best_value = app.m_best_value; + m_best_error = app.m_best_error; m_diff_ulimit = app.m_diff_ulimit; m_flux_ulimit = app.m_flux_ulimit; m_eflux_ulimit = app.m_eflux_ulimit; @@ -603,24 +608,17 @@ void ctulimit::get_parameter_brackets(double& parmin, double& parmax) log_header1(TERSE, "Find parameter brackets"); // Extract current value and error - double value = m_model_par->factor_value(); - double error = m_model_par->factor_error(); - - // Re-compute best log-likelihood for value fixed at optimum - m_best_logL = evaluate(*m_model_par, value); - - // Write best log-likelihood - log_value(NORMAL, "Best log-likelihood", - "logL("+gammalib::str(value)+")="+gammalib::str(m_best_logL)); + m_best_value = m_model_par->factor_value(); + m_best_error = m_model_par->factor_error(); // If parameter error is zero then take parameter value as error - if (error == 0.0) { - error = value; + if (m_best_error == 0.0) { + m_best_error = m_best_value; } // Compute parameter bracketing - parmin = value - m_sigma_min * error; - parmax = value + m_sigma_max * error; + parmin = m_best_value - m_sigma_min * m_best_error; + parmax = m_best_value + m_sigma_max * m_best_error; if (m_model_par->has_min() && m_model_par->factor_min() > parmin) { parmin = m_model_par->factor_min(); } @@ -693,8 +691,9 @@ void ctulimit::ulimit_bisection(const double& parmin, const double& parmax) double wrk_min = parmin; double wrk_max = parmax; - // Initialise iteration counter - int iter = 0; + // Initialise iteration counter and restart flag + int iter = 0; + bool restart = false; // Loop until breaking condition is reached while (true) { @@ -715,29 +714,71 @@ void ctulimit::ulimit_bisection(const double& parmin, const double& parmax) double eval_mid = logL - (m_best_logL + m_dlogL); // Log information - log_value(EXPLICIT, "Iteration "+gammalib::str(iter), + log_value(NORMAL, "Iteration "+gammalib::str(iter), "["+gammalib::str(wrk_min)+", "+gammalib::str(wrk_max)+"] " - "dlogL("+gammalib::str(mid)+")="+gammalib::str(eval_mid)); + "logL("+gammalib::str(mid)+")="+gammalib::str(logL)+" " + "dlogL="+gammalib::str(eval_mid)); // Check for convergence inside tolerance if (std::abs(eval_mid) < m_tol) { break; } + // Assess status + int status = 0; + + // Check whether mid-point is close to best parameter value + double eps = (mid != 0.0) ? (mid - m_best_value)/mid : mid - m_best_value; + if (std::abs(eps) < 1.0e-6) { + std::string msg = + " *** WARNING: Parameter range mid-point "+gammalib::str(mid)+" is " + "close to best factor\n" + " value "+gammalib::str(m_best_value)+" yet the " + "log-likelihood value "+gammalib::str(logL)+"\n" + " at the mid-point differs significantly from the " + "log-likelihood value\n" + " "+gammalib::str(m_best_logL)+" for the best factor " + "value."; + log_string(TERSE, msg); + status = 1; + } + // Check for vanishing parameter range - if (std::abs(wrk_min-wrk_max) < 1.0e-6) { - std::string msg = "*** WARNING: Parameter range " - "["+gammalib::str(wrk_min)+", "+ - gammalib::str(wrk_max)+"] has reduced to a " - "narrow interval without reaching convergence!\n" - " The best log-likelihood value is " - "probably not an absolute but only a local " - "minimum.\n" - " Stop iterations."; + eps = (mid != 0.0) ? (wrk_min - wrk_max)/mid : wrk_min - wrk_max; + if (std::abs(eps) < 1.0e-6) { + std::string msg = + " *** WARNING: Parameter range ["+gammalib::str(wrk_min)+", "+ + gammalib::str(wrk_max)+"] has reduced\n" + " to a narrow interval without reaching convergence! " + "The best\n" + " log-likelihood value is probably not an absolute " + "but only a local\n" + " minimum."; log_string(TERSE, msg); - break; + status = 2; } + // Restart bisection if a restart is required + if (!restart && status != 0) { + + // Adopt current logL as best log-likelihood + restart = true; + iter = 0; + wrk_min = parmin; + wrk_max = parmax; + m_best_logL = logL; + + // Signal restart + std::string msg = + " Adopt "+gammalib::str(logL)+" as best log-likelihood and restart " + "bisection"; + log_string(TERSE, msg); + + // Restart + continue; + + } // endif: restart bisection if required + // Change boundaries for further iteration if (eval_mid > 0.0) { wrk_max = mid; // logL too large, new range = [wrk_min, mid] diff --git a/src/ctulimit/ctulimit.hpp b/src/ctulimit/ctulimit.hpp index 6ec756a..0358f74 100755 --- a/src/ctulimit/ctulimit.hpp +++ b/src/ctulimit/ctulimit.hpp @@ -1,7 +1,7 @@ /*************************************************************************** * ctulimit - Upper limit calculation tool * * ----------------------------------------------------------------------- * - * copyright (C) 2015-2018 by Michael Mayer * + * copyright (C) 2015-2021 by Michael Mayer * * ----------------------------------------------------------------------- * * * * This program is free software: you can redistribute it and/or modify * @@ -98,6 +98,8 @@ protected: GModelSky* m_skymodel; //!< Pointer to sky model GModelPar* m_model_par; //!< Pointer to model parameter double m_best_logL; //!< Best fit log likelihood of given model + double m_best_value; //!< Best parameter value factor + double m_best_error; //!< Best parameter value error double m_diff_ulimit; //!< Differential upper limit value double m_flux_ulimit; //!< Flux upper limit value double m_eflux_ulimit; //!< Energy flux upper limits -- libgit2 0.21.4