Add "refit_if_failed" parameter and set it by default to yes.

The "refit_if_failed" parameter instructs ctlike to perform a refit in case of a fit failure. Fit failures are either a stalled fit, an exhaustion of the number of iterations, or a significant difference between the number of observed and predicted events.

Add "refit_if_failed" parameter and set it by default to yes.
The "refit_if_failed" parameter instructs ctlike to perform a refit in case of a fit failure. Fit failures are either a stalled fit, an exhaustion of the number of iterations, or a significant difference between the number of observed and predicted events.
Jürgen Knödlseder
1 parent 10e1f7cc
Showing 6 changed files with 67 additions and 12 deletions
src/ctlike/ctlike.cpp
src/ctlike/ctlike.hpp
src/ctlike/ctlike.par
src/ctulimit/ctulimit.cpp
test/test_csworkflow.py
test/test_ctulimit.py
@@ -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<GChatter>(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<GChatter>((*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;
+}
@@ -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
@@ -33,6 +33,7 @@ outcovmat, f, h, NONE,,, &quot;Output covariance matrix FITS file&quot;
 #========================
 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?"
@@ -709,11 +709,6 @@ void ctulimit::ulimit_bisection(const double&amp; parmin, const double&amp; 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) {
@@ -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(),
@@ -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')