Merge branch 'upper-limit' into integration

Jürgen Knödlseder
2 parents 6755e796 d223acf4
Showing 1 changed file with 31 additions and 3 deletions
src/ctulimit/ctulimit.cpp
 /***************************************************************************
  *                   ctulimit - Upper limit calculation tool               *
  * ----------------------------------------------------------------------- *
- *  copyright (C) 2015-2020 by Michael Mayer                               *
+ *  copyright (C) 2015-2021 by Michael Mayer                               *
  * ----------------------------------------------------------------------- *
  *                                                                         *
  *  This program is free software: you can redistribute it and/or modify   *
@@ -268,7 +268,8 @@ void ctulimit::run(void)
     } // endif: likelihood was zero
-    // Get parameter brackets
+    // Get parameter brackets. Note that this method recomputes m_best_logL
+    // to make sure that the value is coherent for the upper limit search.
     double parmin;
     double parmax;
     get_parameter_brackets(parmin, parmax);
@@ -280,6 +281,7 @@ void ctulimit::run(void)
     log_value(NORMAL, "Model name", m_skymodel->name());
     log_value(NORMAL, "Parameter name", m_model_par->name());
     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",
               "["+gammalib::str(parmin)+", "+gammalib::str(parmax)+"]");
@@ -604,6 +606,13 @@ void ctulimit::get_parameter_brackets(double&amp; parmin, double&amp; parmax)
     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));
+
     // If parameter error is zero then take parameter value as error
     if (error == 0.0) {
         error = value;
@@ -671,7 +680,12 @@ void ctulimit::get_parameter_brackets(double&amp; parmin, double&amp; parmax)
  *            Exceeded maximum number of iterations in the upper limit
  *            search.
  *
- * Calculates the upper limit using a bisection method.
+ * Calculates the upper limit using a bisection method. The bisection is
+ * stopped if the log-likelihood value is within the tolerance of the target
+ * value. The bisection is also stopped putting a warning in the log file if
+ * the parameter interval becomes smaller than a tolerance of 1.0e-6. The
+ * method stops with an exception if the maximum number of iterations is
+ * exhausted.
  ***************************************************************************/
 void ctulimit::ulimit_bisection(const double& parmin, const double& parmax)
 {
@@ -710,6 +724,20 @@ void ctulimit::ulimit_bisection(const double&amp; parmin, const double&amp; parmax)
             break;
         }
+        // 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.";
+            log_string(TERSE, msg);
+            break;
+        }
+
         // Change boundaries for further iteration
         if (eval_mid > 0.0) {
             wrk_max = mid;     // logL too large, new range = [wrk_min, mid]