cterror implementation

Florent Forest
1 parent 6577d6a2
Showing 11 changed files with 331 additions and 162 deletions
configure.ac
pyext/cterror.i
pyext/ctools/obsutils.py
pyext/ctools/tools.i
scripts/cspull.py
src/Makefile.am
src/cterror/Makefile.am
src/cterror/cterror.cpp
src/cterror/cterror.hpp
src/cterror/cterror.par
src/cterror/main.cpp
@@ -890,6 +890,7 @@ AC_CONFIG_FILES([Makefile
                  src/cttsmap/Makefile
                  src/ctbutterfly/Makefile
                  src/ctulimit/Makefile
+                 src/cterror/Makefile
                  src/support/Makefile
                  src/ctools-setup
                  test/Makefile
+/***************************************************************************
+ *                   cterror - Parameter error calculation tool            *
+ * ----------------------------------------------------------------------- *
+ *  copyright (C) 2015 by Florent Forest                                   *
+ * ----------------------------------------------------------------------- *
+ *                                                                         *
+ *  This program is free software: you can redistribute it and/or modify   *
+ *  it under the terms of the GNU General Public License as published by   *
+ *  the Free Software Foundation, either version 3 of the License, or      *
+ *  (at your option) any later version.                                    *
+ *                                                                         *
+ *  This program is distributed in the hope that it will be useful,        *
+ *  but WITHOUT ANY WARRANTY; without even the implied warranty of         *
+ *  MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the          *
+ *  GNU General Public License for more details.                           *
+ *                                                                         *
+ *  You should have received a copy of the GNU General Public License      *
+ *  along with this program.  If not, see <http://www.gnu.org/licenses/>.  *
+ *                                                                         *
+ ***************************************************************************/
+/**
+ * @file cterror.i
+ * @brief Parameter error calculation tool interface definition
+ * @author Florent Forest
+ */
+
+%{
+/* Put headers and other declarations here that are needed for compilation */
+#include "cterror.hpp"
+%}
+
+
+/***********************************************************************//**
+ * @class cterror
+ *
+ * @brief Parameter error calculation tool
+ *
+ ***************************************************************************/
+class cterror : public ctool {
+
+    public:
+    // Constructors and destructors
+    cterror(void);
+    explicit cterror(const GObservations& obs);
+    cterror(int argc, char *argv[]);
+    cterror(const cterror& app);
+    virtual ~cterror(void);
+
+    // Methods
+    void                 clear(void);
+    void                 run(void);
+    void                 save(void);
+    const GObservations& obs(void) const;
+    const GOptimizer*    opt(void) const;
+
+};
+
+/***********************************************************************//**
+ * @brief Parameter error calculation tool Python extensions
+ ***************************************************************************/
+%extend cterror {
+    cterror copy() {
+        return (*self);
+    }
+}
@@ -158,6 +158,37 @@ def fit(obs, log=False, debug=False, edisp=False):
     # Return observations
     return like
+# ================ #
+# Run cterror      #
+# ================ #
+def cterror(obs, log=False, debug=False):
+    """
+    Perform maximum likelihood fitting of observations in the container.
+    
+    Parameters:
+     obs   - Observation container
+    Keywords:
+     log   - Create log file(s)
+     debug - Create screen dump
+     edisp - Apply energy dispersion?
+    """
+    # Allocate cterror application
+    error = ctools.cterror(obs)
+    
+    # Optionally open the log file
+    if log:
+        error.logFileOpen()
+    
+    # Optionally switch-on debugging model
+    if debug:
+        error["debug"].boolean(True)
+
+    # Run cterror application.
+    error.run()
+    
+    # Return observations
+    return error
+
 # ================= #
 # Create counts map #
@@ -59,3 +59,4 @@
 %include "cttsmap.i"
 %include "ctbutterfly.i"
 %include "ctulimit.i"
+%include "cterror.i"
@@ -297,8 +297,9 @@ class cspull(ctools.cscript):
             self.log("\n")
         # Fit model
-        like = obsutils.fit(obs, edisp=self.m_edisp, \
-                            log=self.m_log, debug=self.m_debug)
+        #like = obsutils.fit(obs, edisp=self.m_edisp, \
+        #                    log=self.m_log, debug=self.m_debug)
+        like  = obsutils.cterror(obs, log=self.m_log, debug=self.m_debug)
         # Store results
         logL   = like.opt().value()
@@ -44,8 +44,8 @@ libctools_la_LIBADD = support/libsupport.la \
                       ctcubemask/libctcubemask.la \
                       cttsmap/libcttsmap.la \
                       ctbutterfly/libctbutterfly.la \
-                      ctulimit/libctulimit.la               
-                      
+                      ctulimit/libctulimit.la \
+                      cterror/libcterror.la
 # Distribute scripts
 EXTRA_DIST   = $(scripts)
@@ -9,10 +9,15 @@ bin_PROGRAMS = cterror
 # Sources
 sources = cterror.cpp
+# Libtool library
+noinst_LTLIBRARIES = libcterror.la
+libcterror_la_SOURCES = $(sources)
+
 # Program definition
 cterror_SOURCES = main.cpp
 cterror_LDFLAGS = $(LIB_LINK_OPTION) @LDFLAGS@
 cterror_CPPFLAGS = @CPPFLAGS@
+cterror_LDADD = libcterror.la $(top_srcdir)/src/support/libsupport.la
 # Parameter files
 pfilesdir = $(prefix)/syspfiles
@@ -28,15 +28,16 @@
 #ifdef HAVE_CONFIG_H
 #include <config.h>
 #endif
-#include <stdio>
+#include <iostream>
 #include "cterror.hpp"
 #include "GTools.hpp"
 #include "GOptimizer.hpp"
 /* __ Method name definitions ____________________________________________ */
 #define G_GET_PARAMETERS                         "cterror::get_parameters()"
-#define G_GET_MODEL_PARAMETER               "cterror::get_model_parameter()"
+#define G_ERR_BISECTION         "cterror::error_bisection(double&, double&)"
 #define G_EVALUATE                              "cterror::evaluate(double&)"
+#define G_RUN                                               "cterror::run()"
 /* __ Debug definitions __________________________________________________ */
@@ -230,92 +231,161 @@ void cterror::run(void)
     // Save original models
     GModels models_orig = m_obs.models();
-    // Save original log-likelihood. If the value is zero it has never been
-    // computed hence we compute it now. 
+    // Allocate optimizer
+    m_opt = new GOptimizerLM();
+    // Optimize and save best log-likelihood
+    m_obs.optimize(*m_opt);
+    m_obs.errors(*m_opt);
     m_best_logL = m_obs.logL();
-    if (m_best_logL == 0.0) {
-
-        // Write header
-        if (logTerse()) {
-            log << std::endl;
-            log.header1("Compute best-fit likelihood");
-        }
-
-        // Reoptimize if likelihood was not given before
-        GOptimizerLM* opt = new GOptimizerLM();
-        m_obs.optimize(*opt);
-        m_obs.errors(*opt);
-        m_best_logL = m_obs.logL();
-
-        // Write optimised model into logger
-        if (logTerse()) {
-            log << m_obs.models() << std::endl;
-        }
-
-    } // endif: likelihood was zero
-
-    // Compute parameter bracketing
-    double value  = m_model_par->value();
-    double error  = m_model_par->error();
-    double parmin = value + m_sigma_min * error;
-    double parmax = value + m_sigma_max * error;
     // Write header
     if (logTerse()) {
         log << std::endl;
-        log.header1("Compute upper limit");
-        log << gammalib::parformat("Model name");
-        log << m_skymodel->name() << std::endl;
-        log << gammalib::parformat("Parameter name");
-        log << m_model_par->name() << std::endl;
-        log << gammalib::parformat("Confidence level");
-        log << m_confidence*100.0 << "%" << std::endl;
-        log << gammalib::parformat("Log-likelihood difference");
-        log << m_dlogL << std::endl;
-        log << gammalib::parformat("Initial parameter range");
-        log << "[";
-        log << parmin;
-        log << ", ";
-        log << parmax;
-        log << "]" << std::endl;
+        log.header1("Compute best-fit likelihood");
     }
-    // Compute lower boundary
-    error_bisection(parmin, value);
-    double value_lo = m_model_par->value();
-
+    // Write optimised model into logger
     if (logTerse()) {
-        log << gammalib::parformat("Final lower parameter value");
-        log << value_lo << std::endl;
+        log << m_obs.models() << std::endl;
     }
-    // Reset parameter to previous value
-    m_model_par->value(value);
-    // Compute upper boundary
-    error_bisection(value, parmax);
-    double value_hi = m_model_par->value();
+    // Continue only if source model exists
+    if (m_obs.models().contains(m_srcname)) {
+
+        GModels& models = const_cast<GModels&>(m_obs.models());
+        m_skymodel      = dynamic_cast<GModelSky*>(models[m_srcname]);
+        if (m_skymodel == NULL) {
+            std::string msg = "Source \""+m_srcname+"\" is not a sky model. "
+                    "Please specify the name of a sky model for "
+                    "parameter error computation.";
+            throw GException::invalid_value(G_RUN, msg);
+        }
+
+        int npars = m_skymodel->spectral()->size();
+
+
+        for (int i = 0; i < npars; ++i)
+        {
+
+            // Re-optimize
+            m_obs.optimize(*m_opt);
+
+            if(m_skymodel->spectral()->at(i).is_fixed())
+                continue;
+
+            m_model_par = &(m_skymodel->spectral()->at(i));
+
+            // Compute parameter bracketing
+            m_value  = m_model_par->factor_value();
+            double parmin = std::max(m_model_par->factor_min(),
+                                     m_value - 10*m_model_par->factor_error());
+            double parmax = std::min(m_model_par->factor_max(),
+                                     m_value + 10*m_model_par->factor_error());
+
+            // Write header
+            if (logTerse()) {
+                log << std::endl;
+                log.header1("Compute parameter error");
+                log << gammalib::parformat("Model name");
+                log << m_skymodel->name() << std::endl;
+                log << gammalib::parformat("Parameter name");
+                log << m_model_par->name() << std::endl;
+                log << gammalib::parformat("Confidence level");
+                log << m_confidence*100.0 << "%" << std::endl;
+                log << gammalib::parformat("Log-likelihood difference");
+                log << m_dlogL << std::endl;
+                log << gammalib::parformat("Initial parameter range");
+                log << "[";
+                log << parmin;
+                log << ", ";
+                log << parmax;
+                log << "]" << std::endl;
+            }
+
+            // Compute lower boundary
+            error_bisection(parmin, m_value);
+            double value_lo = m_model_par->factor_value();
+
+            if (logTerse()) {
+                log << gammalib::parformat("Final lower parameter value");
+                log << value_lo << std::endl;
+            }
+
+            // Reset parameter to original value
+            m_model_par->factor_value(m_value);
+
+            // Compute upper boundary
+            error_bisection(m_value, parmax);
+            double value_hi = m_model_par->factor_value();
+
+            if (logTerse()) {
+                log << gammalib::parformat("Final upper parameter");
+                log << value_hi << std::endl;
+            }
+
+            // Reset parameter to original value
+            m_model_par->factor_value(m_value);
+
+            // Compute errors
+            m_error = 0.5*(value_hi - value_lo);
+            double error_max = std::max(value_hi-m_value, m_value-value_lo);
+            double error_min = std::min(value_hi-m_value, m_value-value_lo);
+
+
+            // Write results to logfile
+            if (logTerse()) {
+                log << std::endl;
+                log.header1("Error results");
+                log << gammalib::parformat("Max error");
+                log << error_max;
+                log << gammalib::parformat("Min error");
+                log << error_min;
+                log << gammalib::parformat("Mean parameter error");
+                log << m_error;
+                log << " " << m_model_par->unit();
+                log << std::endl;
+                log << gammalib::parformat("Parameter scale");
+                log << m_model_par->scale();
+                log << std::endl;
+            }
+
+            // Save error result
+            m_model_par->factor_error(m_error);
+
+        }
-    if (logTerse()) {
-        log << gammalib::parformat("Final upper parameter");
-        log << value_hi << std::endl;
     }
-    // Compute error
-    double error = value_hi - value_lo;
-    // Write results to logfile
+    // Return
+    return;
+}
+
+/***********************************************************************//**
+ * @brief Save maps
+ *
+ * This method saves the error to an ascii file.
+ *
+ * @todo No yet implemented as we have no clear use case yet for saving
+ * the result.
+ ***************************************************************************/
+void cterror::save(void)
+{
+    // Write header
     if (logTerse()) {
         log << std::endl;
-        log.header1("Error results");
-        log << gammalib::parformat("Parameter error");
-        log << error;
-        log << " " << m_model_par->unit();
-        log << std::endl;
+        log.header1("Save results");
     }
-    // Recover original models
-    m_obs.models(models_orig);
+
+    // Get output filename
+    std::string outmodel = (*this)["outfile"].filename();
+
+     // Write results out as XML model
+    if (gammalib::toupper(outmodel) != "NONE") {
+        m_obs.models().save(outmodel);
+    }
     // Return
     return;
@@ -334,11 +404,13 @@ void cterror::init_members(void)
 {
     // Initialise user parameters
     m_srcname.clear();
-    m_confidence = 0.95;
+    m_confidence = 0.68;
     m_sigma_min  = 0.0;
     m_sigma_max  = 0.0;
-    m_tol        = 1.0e-6;
+    m_tol        = 1.0e-3;
     m_max_iter   = 50;
+    m_value      = 0.0;
+    m_error      = 0.0;
     // Initialise protected members
     m_obs.clear();
@@ -367,16 +439,11 @@ void cterror::copy_members(const cterror&amp; app)
     m_tol        = app.m_tol;
     m_max_iter   = app.m_max_iter;
-
     // Copy protected members
     m_obs          = app.m_obs;
     m_dlogL        = app.m_dlogL;
     m_best_logL    = app.m_best_logL;
-
-    // Extract model parameter
-    get_model_parameter();
-
     // Return
     return;
 }
@@ -416,7 +483,7 @@ void cterror::get_parameters(void)
     } // endif: there was no observation in the container
-    // If there is are no models associated with the observations then
+    // If there are no models associated with the observations then
     // load now the model definition
     if (m_obs.models().size() == 0) {
@@ -437,11 +504,9 @@ void cterror::get_parameters(void)
         throw GException::invalid_value(G_GET_PARAMETERS, msg);
     }
-    // Get relevant model and parameter for upper limit computation
-    get_model_parameter();
-
     // Get confidence level and transform into log-likelihood difference
     m_confidence = (*this)["confidence"].real();
+
     double sigma = gammalib::erfinv(m_confidence) * gammalib::sqrt_two;
     m_dlogL      = (sigma*sigma) / 2.0;
@@ -449,70 +514,10 @@ void cterror::get_parameters(void)
     m_sigma_min = (*this)["sigma_min"].real();
     m_sigma_max = (*this)["sigma_max"].real();
-    // Read energy values
-    m_eref = (*this)["eref"].real();
-    m_emin = (*this)["emin"].real();
-    m_emax = (*this)["emax"].real();
-
     // Read precision
     m_tol      = (*this)["tol"].real();
     m_max_iter = (*this)["max_iter"].integer();
-    // Return
-    return;
-}
-
-
-/***********************************************************************//**
- * @brief Get application parameters
- *
- * @exception GException::invalid_value
- *            Did not find a valid model parameter
- *
- * Extracts a pointer to the sky model (m_skymodel) and a pointer to the
- * relevant model parameter (m_model_par) from the model container.
- ***************************************************************************/
-void cterror::get_model_parameter(void)
-{
-    // Continue only if source model exists
-    if (m_obs.models().contains(m_srcname)) {
-
-        // Get relevant model and parameter for upper limit computation.
-        GModels& models = const_cast<GModels&>(m_obs.models());
-        m_skymodel      = dynamic_cast<GModelSky*>(models[m_srcname]);
-        if (m_skymodel == NULL) {
-            std::string msg = "Source \""+m_srcname+"\" is not a sky model. "
-                              "Please specify the name of a sky model for "
-                              "upper limit 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);
-        }
-        if (m_skymodel->spectral()->has_par("Normalization")) {
-            m_model_par = &(m_skymodel->spectral()->operator[]("Normalization"));
-        }
-        else if (m_skymodel->spectral()->has_par("Prefactor")) {
-            m_model_par = &(m_skymodel->spectral()->operator[]("Prefactor"));
-        }
-        else if (m_skymodel->spectral()->has_par("Integral")) {
-            m_model_par = &(m_skymodel->spectral()->operator[]("Integral"));
-        }
-        else if (m_skymodel->spectral()->has_par("Value")) {
-            m_model_par = &(m_skymodel->spectral()->operator[]("Value"));
-        }
-        else {
-            std::string msg = "Require spectral parameter \"Normalization\", "
-                              "\"Prefactor\", \"Integral\" or \"Value\" for "
-                              "upper limit computation. The specified source "
-                              "\""+m_srcname+"\" does not have such a parameter.";
-            throw GException::invalid_value(G_GET_MODEL_PARAMETER, msg);
-        }
-
-    } // endif: source model existed
     // Return
     return;
@@ -551,10 +556,35 @@ void cterror::error_bisection(const double&amp; min, const double&amp; 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)+
-                              " has been reached. You may consider to increase"
-                              " the \"max_iter\" parameter and re-run cterror.";
-            throw GException::invalid_value(G_UL_BISECTION, msg);
+            if(wrk_min - m_model_par->factor_min() < m_tol)
+            {
+                std::string msg = "The "+m_model_par->name()+" parameter minimum has been"
+                                  " reached during error calculation. To obtain accurate "
+                                  " errors, consider to set the minimum to a lower value,"
+                                  " and re-run cterror.";
+                if (logTerse()) {
+                    log << msg;
+                }
+                break;
+            }
+            else if(m_model_par->factor_max() - wrk_max < m_tol)
+            {
+                std::string msg = "The "+m_model_par->name()+" parameter maximum has been"
+                                  " reached during error calculation. To obtain accurate "
+                                  " errors, consider to set the maxmimum to a higher value"
+                                  ", and re-run cterror.";
+                if (logTerse()) {
+                    log << msg;
+                }
+                break;
+            }
+            else
+            {
+                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 cterror.";
+                throw GException::invalid_value(G_ERR_BISECTION, msg);
+            }
         }
         // Compute center of boundary
@@ -568,12 +598,27 @@ void cterror::error_bisection(const double&amp; min, const double&amp; max)
             break;
         }
-        // Change boundaries for further iteration
-        if (eval_mid > 0.0) {
-            wrk_max = mid;
+        // If we are on the crescent side of the parabola
+        if (mid > m_value)
+        {
+            // Change boundaries for further iteration
+            if (eval_mid > 0.0) {
+                wrk_max = mid;
+            }
+            else if (eval_mid < 0.0) {
+                wrk_min = mid;
+            }
         }
-        else if (eval_mid < 0.0) {
-            wrk_min = mid;
+        // If we are on the decrescent side of the parabola
+        else
+        {
+            // Change boundaries for further iteration
+            if (eval_mid > 0.0) {
+                wrk_min = mid;
+            }
+            else if (eval_mid < 0.0) {
+                wrk_max = mid;
+            }
         }
         // Increment counter
@@ -602,10 +647,16 @@ double cterror::evaluate(const double&amp; value)
     double logL = 0.0;
     // Check if given parameter is within boundaries
-    if (value > m_model_par->min() && value < m_model_par->max()) {
+    if (value > m_model_par->factor_min() && value < m_model_par->factor_max()) {
         // Change parameter factor
-        m_model_par->value(value);
+        m_model_par->factor_value(value);
+
+        //Fix parameter
+        m_model_par->fix();
+
+        // Re-optimize
+        m_obs.optimize(*m_opt);
         // Evaluate likelihood for new model container
         m_obs.eval();
@@ -613,6 +664,9 @@ double cterror::evaluate(const double&amp; value)
         // Retrieve likelihood
         logL = m_obs.logL();
+        // Free parameter
+        m_model_par->free();
+
     } // endif: value was inside allowed range
     // ... otherwise signal that the parameter went outside the boundaries
@@ -627,6 +681,7 @@ double cterror::evaluate(const double&amp; value)
     // Compute function value
     double logL_difference = logL - m_best_logL - m_dlogL;
+
     // Return
     return logL_difference;
 }
@@ -61,7 +61,7 @@ public:
     void                 run(void);
     void                 save(void);
     const GObservations& obs(void) const;
-
+    const GOptimizer*    opt(void) const;
 protected:
     // Protected methods
@@ -69,7 +69,6 @@ protected:
     void   copy_members(const cterror& app);
     void   free_members(void);
     void   get_parameters(void);
-    void   get_model_parameter(void);
     double evaluate(const double& value);
     void   error_bisection(const double& min, const double& max);
@@ -80,6 +79,8 @@ protected:
     double        m_sigma_max;    //!< Starting value maximum (multiple fit errors above fit values)
     double        m_tol;          //!< Tolerance for limit determination
     int           m_max_iter;     //!< Maximum number of iterations
+    double        m_value;        //!< Parameter value 
+    double        m_error;        //!< Parameter error
     // Protected members
     GObservations m_obs;          //!< Observation container
@@ -87,6 +88,7 @@ protected:
     GModelSky*    m_skymodel;     //!< Pointer to sky model
     GModelPar*    m_model_par;    //!< Pointer to model parameter
     double        m_best_logL;    //!< Best fit log likelihood of given model
+    GOptimizerLM* m_opt;
 };
@@ -102,5 +104,15 @@ const GObservations&amp; cterror::obs(void) const
     return m_obs;
 }
+/***********************************************************************//**
+ * @brief Return optimizer
+ *
+ * @return Pointer to optimizer
+ ***************************************************************************/
+inline
+const GOptimizer* cterror::opt(void) const
+{
+    return m_opt;
+}
 #endif /* CTERROR_HPP */
@@ -19,6 +19,7 @@
 #===================
 inobs,    f, a, "events.fits",,, "Input event list, counts cube or observation definition file"
 inmodel,  f, a, "$CTOOLS/share/models/crab.xml",,, "Source model"
+outfile,  f, a, "results.xml",,, "Save results"
 srcname,  s, a, "Crab",,, "Source of interest"
 expcube,  f, a, "NONE",,, "Exposure cube file (only needed for stacked analysis)"
 psfcube,  f, a, "NONE",,, "PSF cube file (only needed for stacked analysis)"
@@ -32,10 +33,7 @@ irf,      s, a, &quot;cta_dummy_irf&quot;,,, &quot;Instrument response function&quot;
 confidence, r, h, 0.68,0.0,1.0, "Confidence level"
 sigma_min,  r, h, 0.0,,, "Start minimum value (multiple sigmas above fit value)"
 sigma_max,  r, h, 10.0,,, "Start maximum value (multiple sigmas above fit value)"
-eref,       r, h, 1.0,,, "Reference energy for differential limits (TeV)"
-emin,       r, h, 1.0,,, "Minimum energy for flux limits (TeV)"
-emax,       r, h, 100.0,,, "Maximum energy for flux limits (TeV)"
-tol,        r, h, 1e-5,,, "Computation tolerance"
+tol,        r, h, 1e-3,,, "Computation tolerance"
 max_iter,   i, h, 50,,, "Maximum number of iterations"
 #
@@ -28,7 +28,7 @@
 #ifdef HAVE_CONFIG_H
 #include <config.h>
 #endif
-#include <cstdio>
+#include <iostream>
 #include "cterror.hpp"