Correctly compute total number of fit iterations for XML file, fix a bug in the …

…new refitting logic, and revert back to original reference values for the unit tests

Correctly compute total number of fit iterations for XML file, fix a bug in the …
…new refitting logic, and revert back to original reference values for the unit tests
Jürgen Knödlseder
1 parent 573432d6
Showing 5 changed files with 52 additions and 10 deletions
src/ctlike/ctlike.cpp
src/ctlike/ctlike.hpp
src/ctlike/ctlike.par
test/test_csworkflow.py
test/test_ctulimit.py
@@ -430,6 +430,7 @@ void ctlike::init_members(void)
     m_apply_edisp     = false;
     m_fix_spat_for_ts = false;
     m_chatter         = static_cast<GChatter>(2);
+    m_iter            = 0;
     m_logL            = 0.0;
     m_nobs            = 0.0;
     m_npred           = 0.0;
@@ -459,6 +460,7 @@ void ctlike::copy_members(const ctlike&amp; app)
     m_apply_edisp     = app.m_apply_edisp;
     m_fix_spat_for_ts = app.m_fix_spat_for_ts;
     m_chatter         = app.m_chatter;
+    m_iter            = app.m_iter;
     m_logL            = app.m_logL;
     m_nobs            = app.m_nobs;
     m_npred           = app.m_npred;
@@ -557,6 +559,9 @@ void ctlike::optimize_lm(void)
     log_header1(TERSE, "Maximum likelihood optimisation");
     log.indent(1);
  
+    // Initialise number of iterations
+    m_iter = 0;
+
     // Compute number of fitted parameters
     int nfit = 0;
     for (int i = 0; i < m_obs.models().size(); ++i) {
@@ -578,6 +583,9 @@ void ctlike::optimize_lm(void)
     // Perform LM optimization
     m_obs.optimize(m_opt);
  
+    // Add number of iterations
+    m_iter += m_opt.iter();
+
     // Optionally refit
     if (refit(&m_opt)) {
  
@@ -589,6 +597,9 @@ void ctlike::optimize_lm(void)
         // Optimise again
         m_obs.optimize(m_opt);
  
+        // Add number of iterations
+        m_iter += m_opt.iter();
+
     }
  
     // Optionally show curvature matrix
@@ -701,7 +712,7 @@ GXml ctlike::xml_result(void) const
     result->append(GXmlElement("status", status));
     result->append(GXmlElement("log-likelihood", m_logL));
     result->append(GXmlElement("precision", m_opt.eps()));
-    result->append(GXmlElement("iterations", m_opt.iter()));
+    result->append(GXmlElement("iterations", m_iter));
     result->append(GXmlElement("lambda", m_opt.lambda()));
     result->append(GXmlElement("total_parameters", m_opt.npars()));
     result->append(GXmlElement("fitted_parameters", m_opt.nfree()));
@@ -748,12 +759,48 @@ bool ctlike::refit(const GOptimizer* opt)
         }
  
         // ... otherwise were the number of fit iterations exhausted?
-        else if (opt->iter() >= 1) {
+        else if (opt->iter() >= m_max_iter) {
             refit = true;
             log_value(NORMAL, "Refit requested",
                               "Number of fit iterations exhausted");
         }
  
+        // ... otherwise check if there is a significant difference between
+        // the number of observed and predicted events
+        else {
+            double npred      = m_obs.npred();
+            double difference = m_nobs - npred;
+            if (m_nobs > 0.0) {
+                double fraction   = difference / m_nobs;
+                if (std::abs(fraction) > 1.0e-4) {
+                    refit = true;
+                    log_value(NORMAL, "Refit requested",
+                              "Difference "+
+                              gammalib::str(difference)+
+                              " between number of observed events "+
+                              gammalib::str(m_nobs)+
+                              " and predicted events "+
+                              gammalib::str(npred)+
+                              " is larger than +/- 1.0e-4 times the number of"+
+                              " observed events ("+
+                              gammalib::str(fraction)+").");
+                }
+            }
+            else {
+                if (std::abs(difference) > 10.0) {
+                    refit = true;
+                    log_value(NORMAL, "Refit requested",
+                              "Difference "+
+                              gammalib::str(difference)+
+                              " between number of observed events "+
+                              gammalib::str(m_nobs)+
+                              " and predicted events "+
+                              gammalib::str(npred)+
+                              " is larger than +/- 10.");
+                }
+            }
+        } // endelse: check number of events
+
     }
  
     // Return refit flag
@@ -80,6 +80,7 @@ protected:
     GChatter  m_chatter;         //!< Chattiness
  
     // Members
+    int       m_iter;            //!< Number of iterations
     double    m_logL;            //!< Maximum log likelihood
     double    m_nobs;            //!< Number of observed events
     double    m_npred;           //!< Number of predicted events
@@ -33,7 +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?"
+refit_if_failed, b, h, yes,,, "Refit if fit failed?"
 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?"
@@ -142,8 +142,7 @@ class Test(test):
         prefactor        = 3.63152145731529e-16
         index            = 2.39100016863397
         pre_background   = 1.30910556742873
-        #index_background = 0.252909973473968
-        index_background = 0.252938308459647
+        index_background = 0.252909973473968
  
         # Check fit result values
         self.test_value(models['Crab'][2].value(),
@@ -209,11 +209,6 @@ 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')