From d9f31901af4bbdd60f1faef8b38d84dfe3fcd7fe Mon Sep 17 00:00:00 2001 From: Michael Mayer Date: Mon, 9 Feb 2015 10:26:37 +0100 Subject: [PATCH] implemented first version of ctulimit --- doc/source/reference_manual/ctulimit.rst | 95 +++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++ doc/source/reference_manual/index.rst | 1 + pyext/ctulimit.i | 64 ++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++ src/ctulimit/ctulimit.cpp | 312 +++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++------------------------------------------------- src/ctulimit/ctulimit.hpp | 60 ++++++++++++++++++++++++++++++++++++++++++++++++------------ src/ctulimit/ctulimit.par | 27 +++++++++++++++++++-------- src/ctulimit/main.cpp | 2 +- test/test_ctools.sh | 24 ++++++++++++++++++++++++ test/test_ctulimit.py | 94 ++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++ test/test_python.py | 4 ++++ 10 files changed, 613 insertions(+), 70 deletions(-) create mode 100644 doc/source/reference_manual/ctulimit.rst create mode 100644 pyext/ctulimit.i create mode 100755 test/test_ctulimit.py diff --git a/doc/source/reference_manual/ctulimit.rst b/doc/source/reference_manual/ctulimit.rst new file mode 100644 index 0000000..2592ad1 --- /dev/null +++ b/doc/source/reference_manual/ctulimit.rst @@ -0,0 +1,95 @@ +.. _ctulimit: + +ctulimit +=========== + +Computes upper limit. + + +Synopsis +-------- + +Calculates the upper limit flux of a given source and stores the value in ascii file. +The input model has to be optimised for ctulimit to converge. + +General parameters +------------------ + +``inobs [file]`` + Input event list, counts cube or observation definition XML file. + +``inmodel [file]`` + Model XML file containing the optimised source and background definitions. + +``srcname [string]`` + Name of model component for which upper limit should be computed. + +``expcube = NONE [file]`` + Exposure cube file (only needed for stacked analysis). + +``psfcube = NONE [file]`` + PSF cube file (only needed for stacked analysis). + +``caldb [string]`` + Calibration database. + +``irf [string]`` + Instrumental response function. + +``outfile [file]`` + Output ASCII file name. + +``(cl = 0.95) [real]`` + Confidence Level of upper limit + +``(sigma_min = 0.0) [real]`` + Minimum boundary to start searching for upper limit value. Number of standard deviations above best fit value + +``(sigma_max = 10.0) [real]`` + Maximum boundary to start searching for upper limit value. Number of standard deviations above best fit value + +``(emin = 1.0) [real])`` + Minimum energy of flux upper limit (in TeV). + +``(emax = 100) [real]`` + Maximum energy of flux upper limit (in TeV). + +``(tol = 1e-5) [real]`` + Computation tolerance for minimum + +``(max_iter = 50) [integer]`` + Maximum number of iterations before throwing an exception + + +Standard parameters +------------------- + +``(chatter = 2) [integer]`` + Verbosity of the executable: + chatter = 0: no information will be logged + + chatter = 1: only errors will be logged + + chatter = 2: errors and actions will be logged + + chatter = 3: report about the task execution + + chatter = 4: detailed report about the task execution + +``(clobber = yes) [boolean]`` + Specifies whether an existing output counts cube should be overwritten. + +``(debug = no) [boolean]`` + Enables debug mode. In debug mode the executable will dump any log file output to the console. + +``(mode = ql) [string]`` + Mode of automatic parameters (default is "ql", i.e. "query and learn"). + +``(logfile = ctskymap.log) [string]`` + Name of log file. + + +Related tools +------------- + +None diff --git a/doc/source/reference_manual/index.rst b/doc/source/reference_manual/index.rst index e396ddf..76bcf01 100644 --- a/doc/source/reference_manual/index.rst +++ b/doc/source/reference_manual/index.rst @@ -27,6 +27,7 @@ ctools ctselect --- Selects event data ctskymap --- Generates sky map cttsmap --- Generates Test Statistics map + ctulimit --- Calculates upper limit cscripts diff --git a/pyext/ctulimit.i b/pyext/ctulimit.i new file mode 100644 index 0000000..a65cdf4 --- /dev/null +++ b/pyext/ctulimit.i @@ -0,0 +1,64 @@ +/*************************************************************************** + * ctulimit - ctulimit calculation tool * + * ----------------------------------------------------------------------- * + * copyright (C) 2015 by Michael Mayer * + * ----------------------------------------------------------------------- * + * * + * 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 ctulimit + * @brief upper limit calculation tool + * @author Michael Mayer + */ +%{ +/* Put headers and other declarations here that are needed for compilation */ +#include "ctulimit.hpp" +%} + + +/***********************************************************************//** + * @class ctulimit + * + * @brief upper limit calculation tool + ***************************************************************************/ +class ctulimit : public ctool { +public: + // Constructors and destructors + ctulimit(void); + explicit ctulimit(const GObservations& obs); + ctulimit(int argc, char *argv[]); + ctulimit(const ctulimit& app); + virtual ~ctulimit(void); + + // Methods + void clear(void); + void run(void); + void save(void); + const GObservations& obs(void) const; + const double& diff_ulimit(void) const; + const double& flux_ulimit(void) const; + const double& eflux_ulimit(void) const; +}; + + +/***********************************************************************//** + * @brief upper limit calculation tool Python extensions + ***************************************************************************/ +%extend ctulimit { + ctulimit copy() { + return (*self); + } +} diff --git a/src/ctulimit/ctulimit.cpp b/src/ctulimit/ctulimit.cpp index 6925e1b..1fa49d7 100644 --- a/src/ctulimit/ctulimit.cpp +++ b/src/ctulimit/ctulimit.cpp @@ -1,5 +1,5 @@ /*************************************************************************** - * ctulimitt - upper limit calculation tool * + * ctulimit - upper limit calculation tool * * ----------------------------------------------------------------------- * * copyright (C) 2014 by Michael Mayer * * ----------------------------------------------------------------------- * @@ -31,10 +31,12 @@ #include #include "ctulimit.hpp" #include "GTools.hpp" +#include "GOptimizer.hpp" /* __ Method name definitions ____________________________________________ */ #define G_GET_PARAMETERS "ctulimit::get_parameters()" - +#define G_UL_BISECTION "ctulimit::ul_bisection(const double&, const double&)" +#define G_EVALUATE "ctulimit::evaluate(const double&)" /* __ Debug definitions __________________________________________________ */ /* __ Coding definitions _________________________________________________ */ @@ -225,17 +227,200 @@ void ctulimit::run(void) log << m_obs << std::endl; } + GModels models_orig = m_obs.models(); + + + m_bestloglike = m_obs.logL(); + if (m_bestloglike == 0.0) { + + // Write header + if (logTerse()) { + log << std::endl; + log.header1("Recompute best-fit likelihood"); + } + + // Reoptimize if likelihood was not given before + GOptimizerLM* opt = new GOptimizerLM(); + m_obs.optimize(*opt); + m_obs.errors(*opt); + m_bestloglike = m_obs.logL(); + + } //endif: likelihood was 0.0 + + // Store optimized models and pointer to skymodel of interest + m_models = m_obs.models(); + m_skymodel = dynamic_cast(m_models[m_srcname])->clone(); + + // Get reference energy for differential upper limit + GEnergy eref = GEnergy(m_eref, "TeV"); + + // Create energy range for flux limits + GEnergy emin = GEnergy(m_emin, "TeV"); + GEnergy emax = GEnergy(m_emax, "TeV"); + + // Get value, scale and error of parameter + double value = (*m_skymodel->spectral())[0].value(); + double error = (*m_skymodel->spectral())[0].error(); + double scale = (*m_skymodel->spectral())[0].scale(); + + // Compute starting boundaries + double parmin = value / scale + m_sigma_min * error / scale; + double parmax = value / scale + m_sigma_max * error / scale; + // Write header if (logTerse()) { log << std::endl; log.header1("Compute upper limit"); + log << "Searching for upper limit between "; + log << parmin * scale; + log << " and "; + log << parmax * scale; + log << ": "<spectral()->eval(eref, GTime()); + m_flux_ulimit = m_skymodel->spectral()->flux(emin, emax); + m_eflux_ulimit = m_skymodel->spectral()->eflux(emin, emax); + + // Write results to logfile + if (logTerse()) { + log << std::endl; + log.header1("Upper limit computation finished"); + log << "Upper limits are:" << std::endl; + log << " Differential at "; + log << m_eref <<" TeV: "; + log << m_diff_ulimit << " ph/cm2/s/MeV" << std::endl; + log << " Integral ["; + log << m_emin << " - " << m_emax; + log << "] TeV: "; + log << m_flux_ulimit << " ph/cm2/s" < (*m_skymodel->spectral())[0].factor_min() && value < (*m_skymodel->spectral())[0].factor_max()) { + + // Remove source model + m_models.remove(m_srcname); + + // Change spectral model value + (*m_skymodel->spectral())[0].factor_value(value); + + // Re-append model to container + m_models.append(*m_skymodel); + + // Assign new models to observations + m_obs.models(m_models); + + // Evaluate likelihood for new model container + m_obs.eval(); + + // Retrieve likelihood + LogL = m_obs.logL(); + + } // endif: value was inside allowed range + + else { + // throw exception if value out of range + std::string msg = "Value out of range requested: " + "To omit this error, you could increase the allowed " + "parameter range in the model xml file."; + throw GException::invalid_value(G_EVALUATE, msg); } + // Compute function value + double difflike = LogL - m_bestloglike - m_dloglike; - // Compute upper limit here + // Return + return difflike; +} + +/***********************************************************************//** + * @brief Induces upper limit computation by using a bisection method + * + * This method calculates the upper limit using a bisection method + ***************************************************************************/ +void ctulimit::ulimit_bisection(const double& min, const double& max, const double& scale) +{ + // copy values to working values + double wrk_min = min; + double wrk_max = max; + + // Initialise counter + int iter=0; + + // Loop until breaking condition is reached + while (true) { + + // Log information + if (logExplicit()) { + log << " Iteration "; + log << iter; + log << ": Parameter range reduced to ["; + log << wrk_min * scale; + log <<"; "; + log << wrk_max * scale; + log << "]" << std::endl; + } + + // Throw exception if maximum iterations are reached + if( iter > m_max_iter) { + throw GException::invalid_value(G_UL_BISECTION, "Maximum iterations reached");; + } + + // compute center of boundary + double mid = (wrk_min + wrk_max) / 2.0; + + // Calculate function value + double eval_mid = evaluate(mid); + + // Check for convergence inside tolerance + if (std::abs(eval_mid) < m_tol) { + break; + } + + // change boundaries for further iteration + if (eval_mid > 0.0) { + wrk_max = mid; + } + else if (eval_mid < 0.0) { + wrk_min = mid; + } + + // increment counter + iter++; + + } // endwhile // Return return; + } @@ -249,12 +434,21 @@ void ctulimit::save(void) // Write header if (logTerse()) { log << std::endl; - log.header1("Save TS map"); + log.header1("Save upper limit"); } // Get output filename m_outfile = (*this)["outfile"].filename(); + // Create CSV table with 3 columns + GCsv table(1, 3); + + table.real(0, 0, m_diff_ulimit); + table.real(0, 1, m_flux_ulimit); + table.real(0, 2, m_eflux_ulimit); + + // Save CSV table + table.save(m_outfile, " ", clobber()); // Return return; @@ -273,13 +467,25 @@ void ctulimit::save(void) void ctulimit::init_members(void) { // Initialise members - m_infile.clear(); m_outfile.clear(); - + m_srcname.clear(); + m_bestloglike = 0.0; + m_flux_ulimit = 0.0; + m_diff_ulimit = 0.0; + m_eflux_ulimit = 0.0; + m_dloglike = 0.0; + m_tol = 1e-6; + m_max_iter = 50; + m_eref = 0.0; + m_emin = 0.0; + m_emax = 0.0; + m_sigma_min = 0.0; + m_sigma_max = 0.0; + m_skymodel = NULL; // Initialise protected members m_obs.clear(); - + m_models.clear(); // Return return; @@ -294,15 +500,27 @@ void ctulimit::init_members(void) void ctulimit::copy_members(const ctulimit& app) { // Copy attributes - m_infile = app.m_infile; m_outfile = app.m_outfile; - m_modelfile = app.m_modelfile; m_srcname = app.m_srcname; - m_parname = app.m_parname; - m_algorithm = app.m_algorithm; + m_diff_ulimit = app.m_diff_ulimit; + m_flux_ulimit = app.m_flux_ulimit; + m_eflux_ulimit = app.m_eflux_ulimit; + m_bestloglike = app.m_bestloglike; + m_dloglike = app.m_dloglike; + m_tol = app.m_tol; + m_max_iter = app.m_max_iter; + m_eref = app.m_eref; + m_emin = app.m_emin; + m_emax = app.m_emax; + m_sigma_min = app.m_sigma_min; + m_sigma_max = app.m_sigma_max; // Copy protected members m_obs = app.m_obs; + m_models = app.m_models; + + // Clone protected members + m_skymodel = (app.m_skymodel != NULL) ? app.m_skymodel->clone() : NULL; // Return return; @@ -334,51 +552,25 @@ void ctulimit::free_members(void) ***************************************************************************/ void ctulimit::get_parameters(void) { - // If there are no observations in container then add a single CTA - // observation using the parameters from the parameter file + // If there are no observations in container then load them via user + // parameters if (m_obs.size() == 0) { - // Allocate CTA observation - GCTAObservation obs; - - // Get event file name - std::string filename = (*this)["infile"].filename(); - - // Try first to open as FITS file - try { - - // Load data - obs.load(filename); - - // Set response - set_obs_response(&obs); - - // Append observation to container - m_obs.append(obs); - - } - - // ... otherwise try to open as XML file - catch (GException::fits_open_error &e) { - - // Load observations from XML file - m_obs.load(filename); - - // Check if all observations have response information. If - // not, get the calibration database parameters and set - // the response properly - set_response(m_obs); + // Throw exception if no input observation file is given + require_inobs(G_GET_PARAMETERS); - } // endcatch: file was an XML file + // Build observation container + m_obs = get_observations(); } // endif: there was no observation in the container + // If there is are no models associated with the observations then // load now the model definition if (m_obs.models().size() == 0) { // Get models XML filename - std::string filename = (*this)["srcmdl"].filename(); + std::string filename = (*this)["inmodel"].filename(); // Setup models for optimizing. m_obs.models(GModels(filename)); @@ -390,10 +582,34 @@ void ctulimit::get_parameters(void) if (!m_obs.models().contains(m_srcname)) { std::string msg = "Source \""+m_srcname+"\" not found in model " "container. Please add a source with that name " - "or check for a possible typos."; - throw GException::invalid_value(G_GET_PARAMETERS, msg); + "or check for possible typos."; + throw GException::invalid_value(G_GET_PARAMETERS, msg); + } + + double CL = (*this)["cl"].real(); + if (CL != 0.95) { + std::string msg = "Confidence level different from 95% requested" + "Currently only 95% CL is possible"; + throw GException::invalid_value(G_GET_PARAMETERS, msg); + } + else { + // Set Likelihood difference for 95% CL. + // See Minuit Handbook + m_dloglike = 3.84 / 2.0; } + // Read starting boundaries for bisection + 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(); // Optionally read ahead parameters so that they get correctly // dumped into the log file @@ -401,8 +617,6 @@ void ctulimit::get_parameters(void) m_outfile = (*this)["outfile"].filename(); } - // TODO read other parameters here - // Return return; } diff --git a/src/ctulimit/ctulimit.hpp b/src/ctulimit/ctulimit.hpp index b100ea7..f3d0cad 100644 --- a/src/ctulimit/ctulimit.hpp +++ b/src/ctulimit/ctulimit.hpp @@ -69,7 +69,9 @@ public: void run(void); void save(void); const GObservations& obs(void) const; - const double& ulimit(void) const; + const double& diff_ulimit(void) const; + const double& flux_ulimit(void) const; + const double& eflux_ulimit(void) const; protected: // Protected methods @@ -77,19 +79,32 @@ protected: void copy_members(const ctulimit& app); void free_members(void); void get_parameters(void); + void ulimit_bisection(const double& parmin, const double& parmax, const double& scale); + double evaluate(const double& value); // User parameters - std::string m_infile; //!< Input file - std::string m_modelfile;//!< Name of model optimised model file std::string m_srcname; //!< Name of source which is moved around - std::string m_parname; //!< Name of parameter upper limit should be computed - std::string m_outfile; //!< Output counts map or XML file - std::string m_algorithm; //!< algorithm of upper limit computation + std::string m_outfile; //!< Output ascii file + + GModelSky* m_skymodel; //!< Pointer to variable spectral model + + int m_max_iter; //!< maximum number of iterations + double m_dloglike; //!< Likelihood difference for upper limit computation + double m_bestloglike ; //!< Best fit log likelihood of given model + double m_tol; //!< tolerance for limit determination + 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) + double m_eref; //!< reference energy for flux limits + double m_emin; //!< minimum energy for flux limits + double m_emax; //!< maximum energy for flux limits + + double m_diff_ulimit; //!. +# +# ========================================================================== +import gammalib +import ctools + + +# ========================== # +# Test class for cttsmap tool # +# ========================== # +class Test(gammalib.GPythonTestSuite): + """ + Test class for cttsmap tool. + """ + # Constructor + def __init__(self): + """ + Constructor. + """ + # Call base class constructor + gammalib.GPythonTestSuite.__init__(self) + + # Set members + self.events_name = "data/crab_events.fits" + self.model_name = "result.xml" + self.caldb = "irf" + self.irf = "cta_dummy_irf" + + # Return + return + + # Set test functions + def set(self): + """ + Set all test functions. + """ + # Set test name + self.name("ctulimit") + + # Append tests + self.append(self.test_functional, "Test ctulimit functionality") + + # Return + return + + # Test cttsmap functionnality + def test_functional(self): + """ + Test cttsmap functionnality. + """ + # Set-up cttsmap + ulimit = ctools.ctulimit() + ulimit["inobs"].filename(self.events_name) + ulimit["inmodel"].filename(self.model_name) + ulimit["srcname"].string("Crab") + ulimit["caldb"].string(self.caldb) + ulimit["irf"].string(self.irf) + ulimit["outfile"].filename("ulimit.dat") + + # Run tool + self.test_try("Run ctulimit") + try: + ulimit.run() + self.test_try_success() + except: + self.test_try_failure("Exception occured in ctulimit.") + + # Save results + self.test_try("Save results") + try: + ulimit.save() + self.test_try_success() + except: + self.test_try_failure("Exception occured in saving results.") + + # Return + return diff --git a/test/test_python.py b/test/test_python.py index 6dd5aaa..7da0400 100755 --- a/test/test_python.py +++ b/test/test_python.py @@ -34,6 +34,7 @@ import test_ctpsfcube import test_ctbkgcube import test_ctcubemask import test_ctbutterfly +import test_ctulimit import test_pipelines @@ -60,6 +61,7 @@ if __name__ == '__main__': suite_ctbkgcube = test_ctbkgcube.Test() suite_ctcubemask = test_ctcubemask.Test() suite_ctbutterfly = test_ctbutterfly.Test() + suite_ctulimit = test_ctulimit.Test() suite_pipelines = test_pipelines.Test() # Setup unit tests @@ -75,6 +77,7 @@ if __name__ == '__main__': suite_ctbkgcube.set() suite_ctcubemask.set() suite_ctbutterfly.set() + suite_ctulimit.set() suite_pipelines.set() # Append tests to container @@ -90,6 +93,7 @@ if __name__ == '__main__': suites.append(suite_ctbkgcube) suites.append(suite_ctcubemask) suites.append(suite_ctbutterfly) + suites.append(suite_ctulimit) suites.append(suite_pipelines) # Set PFILES environment variable -- libgit2 0.21.4