Add "parname" parameter to ctulimit tool (#2233)

Jürgen Knödlseder
1 parent 5aa6adb3
Showing 4 changed files with 156 additions and 44 deletions
NEWS
src/ctulimit/ctulimit.cpp
src/ctulimit/ctulimit.hpp
src/ctulimit/ctulimit.par
 New Features and Important Changes in ctools 1.6.0
  
-15 June 2018
+20 June 2018
  
  
 Introduction
@@ -151,7 +151,9 @@ None
  
 ctulimit - Compute upper limit
 ------------------------------
-None
+Add "parname" parameter to allow upper limit computations also for a node
+function spectral model. In addition, it is now checked that the parameter
+values actually bracket the upper limit (#2233).
  
  
 support - Support classes
@@ -35,6 +35,8 @@
  
 /* __ Method name definitions ____________________________________________ */
 #define G_GET_MODEL_PARAMETER               "ctulimit::get_model_parameter()"
+#define G_GET_PARAMETER_BRACKETS "ctulimit::get_parameter_brackets(double&, "\
+                                                                   "double&)"
 #define G_UL_BISECTION             "ctulimit::ul_bisection(double&, double&)"
  
 /* __ Debug definitions __________________________________________________ */
@@ -266,24 +268,10 @@ void ctulimit::run(void)
  
     } // endif: likelihood was zero
  
-    // Extract current value and error
-    double value = m_model_par->factor_value();
-    double error = m_model_par->factor_error();
-
-    // If parameter error is zero then take parameter value as error
-    if (error == 0) {
-        error = value;
-    }
-
-    // Compute parameter bracketing
-    double parmin = value - m_sigma_min * error;
-    double parmax = value + m_sigma_max * error;
-    if (m_model_par->has_min() && m_model_par->factor_min() > parmin) {
-        parmin = m_model_par->factor_min();
-    }
-    if (m_model_par->has_max() && m_model_par->factor_max() < parmax) {
-        parmax = m_model_par->factor_max();
-    }
+    // Get parameter brackets
+    double parmin;
+    double parmax;
+    get_parameter_brackets(parmin, parmax);
  
     // Write header into logger
     log_header1(TERSE, "Compute upper limit");
@@ -363,6 +351,7 @@ void ctulimit::init_members(void)
 {
     // Initialise user parameters
     m_srcname.clear();
+    m_parname.clear();
     m_confidence  = 0.95;
     m_sigma_min   = 0.0;
     m_sigma_max   = 0.0;
@@ -401,6 +390,7 @@ void ctulimit::copy_members(const ctulimit&amp; app)
 {
     // Copy user parameters
     m_srcname     = app.m_srcname;
+    m_parname     = app.m_parname;
     m_confidence  = app.m_confidence;
     m_sigma_min   = app.m_sigma_min;
     m_sigma_max   = app.m_sigma_max;
@@ -453,8 +443,9 @@ void ctulimit::get_parameters(void)
     // Setup models from "inmodel" parameter
     setup_models(m_obs, (*this)["srcname"].string());
  
-    // Get name of test source
+    // Get name of test source and optional parameter
     m_srcname = (*this)["srcname"].string();
+    m_parname = (*this)["parname"].string();
  
     // Get relevant model and parameter for upper limit computation
     get_model_parameter();
@@ -532,22 +523,47 @@ void ctulimit::get_model_parameter(void)
                           "computation.";
         throw GException::invalid_value(G_GET_MODEL_PARAMETER, msg);
     }
-    if (m_skymodel->spectral()->type() == "NodeFunction") {
-        std::string msg = "\"NodeFunction\" cannot be used as spectral model "
-                          "for an upper limit computation. Please specify "
-                          "another spectral model.";
-        throw GException::invalid_value(G_GET_MODEL_PARAMETER, msg);
-    }
  
-    // Find appropriate model parameter and store its pointer in the
-    // m_model_par member
-    for (const char** par = pars; *par != NULL; ++par) {
-        if (m_skymodel->spectral()->has_par(*par)) {
-            m_model_par = &(m_skymodel->spectral()->operator[](*par));
-            break;
+    // If a parameter name was specified then use that parameter
+    if (!m_parname.empty()) {
+        if (m_skymodel->spectral()->has_par(m_parname)) {
+            m_model_par = &(m_skymodel->spectral()->operator[](m_parname));
+        }
+        else {
+            std::string msg = "Spectral model of source \""+m_srcname+"\" has "
+                              "no parameter \""+m_parname+"\". Please specify "
+                              "a value parameter name of leave the parameter "
+                              "name blank for autodetermination of an "
+                              "intensity- or flux-like parameter.";
+            throw GException::invalid_value(G_GET_MODEL_PARAMETER, msg);
         }
     }
  
+    // ... otherwise find the appropriate model parameter
+    else {
+
+        // If a node function was specified then a valid parameter name is
+        // needed
+        if (m_skymodel->spectral()->type() == "NodeFunction") {
+            std::string msg = "Spectral model of source \""+m_srcname+"\" is "
+                              "a \"NodeFunction\" but no parameter name was "
+                              "specified. Please use the \"parname\" parameter "
+                              "to specify the parameter that should be used "
+                              "for the upper limit computation.";
+            throw GException::invalid_value(G_GET_MODEL_PARAMETER, msg);
+        }
+
+        // Find appropriate model parameter and store its pointer in the
+        // m_model_par member
+        for (const char** par = pars; *par != NULL; ++par) {
+            if (m_skymodel->spectral()->has_par(*par)) {
+                m_model_par = &(m_skymodel->spectral()->operator[](*par));
+                break;
+            }
+        }
+
+    } // endelse: searched for parameter
+
     // If no model parameter was found then throw an exception
     if (m_model_par == NULL) {
         std::string msg = "Require one of the following spectral "
@@ -569,18 +585,105 @@ void ctulimit::get_model_parameter(void)
  
  
 /***********************************************************************//**
+ * @brief Determine parameter brackets
+ *
+ * @param[out] parmin Minimum parameter value.
+ * @param[out] parmax Maximum parameter value.
+ *
+ * @exception GException::invalid_value
+ *            Exceeded maximum number of iterations in the bracket search.
+ *
+ * Determines the brackets for the parameter of interest.
+ ***************************************************************************/
+void ctulimit::get_parameter_brackets(double& parmin, double& parmax)
+{
+    // Write header into logger
+    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();
+
+    // If parameter error is zero then take parameter value as error
+    if (error == 0) {
+        error = value;
+    }
+
+    // If parameter error is larger than parameter value then take parameter
+    // value as error
+    else if (std::abs(error) > std::abs(value)) {
+        error = value;
+    }
+
+    // Compute parameter bracketing
+    parmin = value - m_sigma_min * error;
+    parmax = value + m_sigma_max * error;
+    if (m_model_par->has_min() && m_model_par->factor_min() > parmin) {
+        parmin = m_model_par->factor_min();
+    }
+    if (m_model_par->has_max() && m_model_par->factor_max() < parmax) {
+        parmax = m_model_par->factor_max();
+    }
+
+    // Initialise iteration counter
+    int iter = 0;
+
+    // Make sure that upper limit is bracketed
+    while (true) {
+
+        // Throw exception if maximum iterations are reached
+        if (iter > m_max_iter) {
+            std::string msg = "The maximum number of "+gammalib::str(m_max_iter)+
+                              " has been reached. You may consider to increase"
+                              " the \"max_iter\" parameter and re-run ctulimit.";
+            throw GException::invalid_value(G_GET_PARAMETER_BRACKETS, msg);
+        }
+
+        // Evaluate logL at upper boundary
+        double logL     = evaluate(*m_model_par, parmax);
+        double eval_mid = logL - (m_best_logL + m_dlogL);
+
+        // Write current parameter range into logger
+        log_value(NORMAL, "Parameter range",
+                  "["+gammalib::str(parmin)+", "+gammalib::str(parmax)+"] "
+                  "logL("+gammalib::str(parmax)+")="+gammalib::str(logL));
+
+        // If the logL increase is not enough, increase maximum parameter
+        // value, otherwise break
+        if (eval_mid < 0.0) {
+            parmax *= 10.0;
+        }
+        else {
+            break;
+        }
+
+        // Increase number of iterations
+        iter++;
+
+    } // endwhile
+
+    // Return
+    return;
+}
+
+
+/***********************************************************************//**
  * @brief Calculate upper limit using a bisection method
  *
- * @param[in] min Minimum parameter value
- * @param[in] max Maximum parameter value
+ * @param[in] parmin Minimum parameter value
+ * @param[in] parmax Maximum parameter value
+ *
+ * @exception GException::invalid_value
+ *            Exceeded maximum number of iterations in the upper limit
+ *            search.
  *
  * Calculates the upper limit using a bisection method.
  ***************************************************************************/
-void ctulimit::ulimit_bisection(const double& min, const double& max)
+void ctulimit::ulimit_bisection(const double& parmin, const double& parmax)
 {
     // Copy values to working values
-    double wrk_min = min;
-    double wrk_max = max;
+    double wrk_min = parmin;
+    double wrk_max = parmax;
  
     // Initialise iteration counter
     int iter = 0;
@@ -588,10 +691,6 @@ void ctulimit::ulimit_bisection(const double&amp; min, const double&amp; max)
     // Loop until breaking condition is reached
     while (true) {
  
-        // Log information
-        log_value(EXPLICIT, "Iteration "+gammalib::str(iter),
-                  "["+gammalib::str(wrk_min)+", "+gammalib::str(wrk_max)+"]");
-
         // Throw exception if maximum iterations are reached
         if (iter > m_max_iter) {
             std::string msg = "The maximum number of "+gammalib::str(m_max_iter)+
@@ -604,15 +703,23 @@ void ctulimit::ulimit_bisection(const double&amp; min, const double&amp; max)
         double mid = (wrk_min + wrk_max) / 2.0;
  
         // Calculate function value
-        double eval_mid = evaluate(*m_model_par, mid) - (m_best_logL + m_dlogL);
+        double logL     = evaluate(*m_model_par, mid);
+        double eval_mid = logL - (m_best_logL + m_dlogL);
+
+        // Log information
+        log_value(EXPLICIT, "Iteration "+gammalib::str(iter),
+                  "["+gammalib::str(wrk_min)+", "+gammalib::str(wrk_max)+"] "
+                  "dlogL("+gammalib::str(mid)+")="+gammalib::str(eval_mid));
  
         // Check for convergence inside tolerance
         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) {
@@ -75,11 +75,13 @@ protected:
     void free_members(void);
     void get_parameters(void);
     void get_model_parameter(void);
+    void get_parameter_brackets(double& parmin, double& parmax);
     void ulimit_bisection(const double& parmin, const double& parmax);
     void compute_ulimit(void);
  
     // User parameters
-    std::string   m_srcname;      //!< Name of source which is moved around
+    std::string   m_srcname;      //!< Name of source for upper limit computation
+    std::string   m_parname;      //!< Name of parameter for upper limit computation
     double        m_confidence;   //!< Confidence level
     double        m_sigma_min;    //!< Starting value minimum (multiple fit errors above fit values)
     double        m_sigma_max;    //!< Starting value maximum (multiple fit errors above fit values)
@@ -19,6 +19,7 @@
 inobs,    f, a, events.fits,,, "Input event list, counts cube or observation definition XML file"
 inmodel,  f, a, $CTOOLS/share/models/crab.xml,,, "Input model definition XML file"
 srcname,  s, a, Crab,,, "Source of interest"
+parname,  s, h, ,,, "Parameter name (optional)"
 expcube,  f, a, NONE,,, "Input exposure cube file"
 psfcube,  f, a, NONE,,, "Input PSF cube file"
 edispcube,f, a, NONE,,, "Input energy dispersion cube file"