From 8ded0f5f0e46fecf704d04be8b839991b68dfbba Mon Sep 17 00:00:00 2001 From: Florent Forest Date: Fri, 19 Jun 2015 14:25:35 +0200 Subject: [PATCH] cterror implementation --- configure.ac | 1 + pyext/cterror.i | 65 +++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++ pyext/ctools/obsutils.py | 31 +++++++++++++++++++++++++++++++ pyext/ctools/tools.i | 1 + scripts/cspull.py | 5 +++-- src/Makefile.am | 4 ++-- src/cterror/Makefile.am | 5 +++++ src/cterror/cterror.cpp | 357 ++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++------------------------------------------------------------------------------------------------------------------------------------------------------- src/cterror/cterror.hpp | 16 ++++++++++++++-- src/cterror/cterror.par | 6 ++---- src/cterror/main.cpp | 2 +- 11 files changed, 331 insertions(+), 162 deletions(-) create mode 100644 pyext/cterror.i diff --git a/configure.ac b/configure.ac index d884d73..9442b21 100644 --- a/configure.ac +++ b/configure.ac @@ -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 diff --git a/pyext/cterror.i b/pyext/cterror.i new file mode 100644 index 0000000..07b9e0d --- /dev/null +++ b/pyext/cterror.i @@ -0,0 +1,65 @@ +/*************************************************************************** + * 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 . * + * * + ***************************************************************************/ +/** + * @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); + } +} diff --git a/pyext/ctools/obsutils.py b/pyext/ctools/obsutils.py index 2ad0d87..4e2181c 100755 --- a/pyext/ctools/obsutils.py +++ b/pyext/ctools/obsutils.py @@ -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 # diff --git a/pyext/ctools/tools.i b/pyext/ctools/tools.i index 381daa8..bc3a0ce 100755 --- a/pyext/ctools/tools.i +++ b/pyext/ctools/tools.i @@ -59,3 +59,4 @@ %include "cttsmap.i" %include "ctbutterfly.i" %include "ctulimit.i" +%include "cterror.i" diff --git a/scripts/cspull.py b/scripts/cspull.py index 9abf3ae..d8aa036 100755 --- a/scripts/cspull.py +++ b/scripts/cspull.py @@ -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() diff --git a/src/Makefile.am b/src/Makefile.am index 2a04273..0709151 100644 --- a/src/Makefile.am +++ b/src/Makefile.am @@ -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) diff --git a/src/cterror/Makefile.am b/src/cterror/Makefile.am index 8ebace8..98cfe05 100644 --- a/src/cterror/Makefile.am +++ b/src/cterror/Makefile.am @@ -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 diff --git a/src/cterror/cterror.cpp b/src/cterror/cterror.cpp index 7e2cbea..95aaf02 100644 --- a/src/cterror/cterror.cpp +++ b/src/cterror/cterror.cpp @@ -28,15 +28,16 @@ #ifdef HAVE_CONFIG_H #include #endif -#include +#include #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(m_obs.models()); + m_skymodel = dynamic_cast(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& 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(m_obs.models()); - m_skymodel = dynamic_cast(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& min, const double& 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& min, const double& 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& 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& 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& value) // Compute function value double logL_difference = logL - m_best_logL - m_dlogL; + // Return return logL_difference; } diff --git a/src/cterror/cterror.hpp b/src/cterror/cterror.hpp index 6c1e0fb..96b3e00 100644 --- a/src/cterror/cterror.hpp +++ b/src/cterror/cterror.hpp @@ -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& 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 */ diff --git a/src/cterror/cterror.par b/src/cterror/cterror.par index 1d5885f..f719705 100644 --- a/src/cterror/cterror.par +++ b/src/cterror/cterror.par @@ -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, "cta_dummy_irf",,, "Instrument response function" 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" # diff --git a/src/cterror/main.cpp b/src/cterror/main.cpp index 895a232..5585a65 100644 --- a/src/cterror/main.cpp +++ b/src/cterror/main.cpp @@ -28,7 +28,7 @@ #ifdef HAVE_CONFIG_H #include #endif -#include +#include #include "cterror.hpp" -- libgit2 0.21.4