From 573432d67780b9c8e56bdafb1b55ab9f07fe6a7f Mon Sep 17 00:00:00 2001 From: Jürgen Knödlseder Date: Fri, 27 Nov 2020 09:45:33 +0100 Subject: [PATCH] Add "refit_if_failed" parameter and set it by default to yes. --- src/ctlike/ctlike.cpp | 59 ++++++++++++++++++++++++++++++++++++++++++++++++++++++----- src/ctlike/ctlike.hpp | 6 +++++- src/ctlike/ctlike.par | 1 + src/ctulimit/ctulimit.cpp | 5 ----- test/test_csworkflow.py | 3 ++- test/test_ctulimit.py | 5 +++++ 6 files changed, 67 insertions(+), 12 deletions(-) diff --git a/src/ctlike/ctlike.cpp b/src/ctlike/ctlike.cpp index 8d88b3b..cf942f6 100755 --- a/src/ctlike/ctlike.cpp +++ b/src/ctlike/ctlike.cpp @@ -423,7 +423,10 @@ void ctlike::init_members(void) // Initialise members m_outmodel.clear(); m_outcovmat.clear(); + m_max_iter = 50; + m_like_accuracy = 0.005; m_refit = false; + m_refit_if_failed = true; m_apply_edisp = false; m_fix_spat_for_ts = false; m_chatter = static_cast(2); @@ -447,9 +450,12 @@ void ctlike::init_members(void) void ctlike::copy_members(const ctlike& app) { // Copy attributes - m_refit = app.m_refit; m_outmodel = app.m_outmodel; m_outcovmat = app.m_outcovmat; + m_max_iter = app.m_max_iter; + m_like_accuracy = app.m_like_accuracy; + m_refit = app.m_refit; + m_refit_if_failed = app.m_refit_if_failed; m_apply_edisp = app.m_apply_edisp; m_fix_spat_for_ts = app.m_fix_spat_for_ts; m_chatter = app.m_chatter; @@ -498,11 +504,14 @@ void ctlike::get_parameters(void) } // endif: no models were associated with observations // Set optimizer characteristics from user parameters - m_opt.eps((*this)["like_accuracy"].real()); - m_opt.max_iter((*this)["max_iter"].integer()); + m_max_iter = (*this)["max_iter"].integer(); + m_like_accuracy = (*this)["like_accuracy"].real(); + m_opt.eps(m_like_accuracy); + m_opt.max_iter(m_max_iter); // Get other parameters m_refit = (*this)["refit"].boolean(); + m_refit_if_failed = (*this)["refit_if_failed"].boolean(); m_apply_edisp = (*this)["edisp"].boolean(); m_fix_spat_for_ts = (*this)["fix_spat_for_ts"].boolean(); m_chatter = static_cast((*this)["chatter"].integer()); @@ -570,7 +579,7 @@ void ctlike::optimize_lm(void) m_obs.optimize(m_opt); // Optionally refit - if (m_refit) { + if (refit(&m_opt)) { // Dump new header log.indent(0); @@ -621,7 +630,7 @@ double ctlike::reoptimize_lm(void) m_obs.optimize(*opt); // Optionally refit - if (m_refit) { + if (refit(opt)) { m_obs.optimize(*opt); } @@ -710,3 +719,43 @@ GXml ctlike::xml_result(void) const // Return XML result return xml; } + + +/***********************************************************************//** + * @brief Refit needed? + * + * @param[in] opt Pointer to optimiser. + * @return True if refit is needed, false otherwise. + * + * This method returns true if either the "refit" parameter is set to yes + * or the "refit_if_failed" parameter is set to yes and the fit has failed. + * The fit is considered as failed if it either has stalled or the number + * of fit iterations is exhausted. + ***************************************************************************/ +bool ctlike::refit(const GOptimizer* opt) +{ + // Initialise refit flag + bool refit = m_refit; + + // If flag is false and "refit_if_failed" is true then examine optimiser + // result to determine whether the fit has failed + if (!refit && m_refit_if_failed) { + + // Has the fit stalled? + if (opt->status() == G_LM_STALLED) { + refit = true; + log_value(NORMAL, "Refit requested", "Fit has stalled"); + } + + // ... otherwise were the number of fit iterations exhausted? + else if (opt->iter() >= 1) { + refit = true; + log_value(NORMAL, "Refit requested", + "Number of fit iterations exhausted"); + } + + } + + // Return refit flag + return refit; +} diff --git a/src/ctlike/ctlike.hpp b/src/ctlike/ctlike.hpp index d3ef320..81917f0 100755 --- a/src/ctlike/ctlike.hpp +++ b/src/ctlike/ctlike.hpp @@ -66,11 +66,15 @@ protected: void optimize_lm(void); double reoptimize_lm(void); GXml xml_result(void) const; + bool refit(const GOptimizer* opt); // User parameters GFilename m_outmodel; //!< Source model output XML file name GFilename m_outcovmat; //!< Covariance matrix output file name - bool m_refit; //!< Refitting + int m_max_iter; //!< Maximum number of iterations + double m_like_accuracy; //!< Log-likelihood accuracy + bool m_refit; //!< Refitting? + bool m_refit_if_failed; //!< Refitting in case of failure? bool m_apply_edisp; //!< Apply energy dispersion? bool m_fix_spat_for_ts; //!< Fix spatial parameters for TS computation? GChatter m_chatter; //!< Chattiness diff --git a/src/ctlike/ctlike.par b/src/ctlike/ctlike.par index 79c5121..00bee25 100755 --- a/src/ctlike/ctlike.par +++ b/src/ctlike/ctlike.par @@ -33,6 +33,7 @@ outcovmat, f, h, NONE,,, "Output covariance matrix FITS file" #======================== statistic, s, h, DEFAULT,DEFAULT|CSTAT|WSTAT|CHI2,, "Optimization statistic" refit, b, h, no,,, "Do refitting?" +refit_if_failed, b, h, yes,,, "Refit if fit stalled or fit iterations were exhausted?" like_accuracy, r, h, 0.005,1.0e-30,1.0, "Absolute accuracy of maximum likelihood value" max_iter, i, h, 50,1,1000, "Maximum number of fit iterations" fix_spat_for_ts, b, h, no,,, "Fix spatial parameters for TS computation?" diff --git a/src/ctulimit/ctulimit.cpp b/src/ctulimit/ctulimit.cpp index b990a99..62df69d 100755 --- a/src/ctulimit/ctulimit.cpp +++ b/src/ctulimit/ctulimit.cpp @@ -709,11 +709,6 @@ void ctulimit::ulimit_bisection(const double& parmin, const double& parmax) if (std::abs(eval_mid) < m_tol) { break; } - /* - if (std::abs(wrk_max-wrk_min) < m_tol) { - break; - } - */ // Change boundaries for further iteration if (eval_mid > 0.0) { diff --git a/test/test_csworkflow.py b/test/test_csworkflow.py index d6cb141..4696432 100755 --- a/test/test_csworkflow.py +++ b/test/test_csworkflow.py @@ -142,7 +142,8 @@ class Test(test): prefactor = 3.63152145731529e-16 index = 2.39100016863397 pre_background = 1.30910556742873 - index_background = 0.252909973473968 + #index_background = 0.252909973473968 + index_background = 0.252938308459647 # Check fit result values self.test_value(models['Crab'][2].value(), diff --git a/test/test_ctulimit.py b/test/test_ctulimit.py index 73fa9a9..4106beb 100755 --- a/test/test_ctulimit.py +++ b/test/test_ctulimit.py @@ -209,6 +209,11 @@ class Test(test): ulimit.logFileOpen() # Needed to get a new log file ulimit.execute() + # Set new reference values since the ctlike results changed + # a bit due to changes in ctlike + ref_flux = 1.6694433147174e-11 + ref_eflux = 6.458852107887e-11 + # Check results self.test_value(ulimit.diff_ulimit(), ref_diff, 1.0e-21, 'Check differential upper limit') -- libgit2 0.21.4