Commit d7ca4cffabd94543ffa44997cc94c327bf5823db
1 parent
bf3c48ec
cterror #1
Showing
5 changed files
with
884 additions
and
0 deletions
src/cterror/Makefile.am
0 → 100644
| 1 | +## Process this file with automake to produce Makefile.in | ||
| 2 | + | ||
| 3 | +# Set include directories | ||
| 4 | +INCLUDES = -I$(top_srcdir)/src/support | ||
| 5 | + | ||
| 6 | +# Program name | ||
| 7 | +bin_PROGRAMS = cterror | ||
| 8 | + | ||
| 9 | +# Sources | ||
| 10 | +sources = cterror.cpp | ||
| 11 | + | ||
| 12 | +# Program definition | ||
| 13 | +cterror_SOURCES = main.cpp | ||
| 14 | +cterror_LDFLAGS = $(LIB_LINK_OPTION) @LDFLAGS@ | ||
| 15 | +cterror_CPPFLAGS = @CPPFLAGS@ | ||
| 16 | + | ||
| 17 | +# Parameter files | ||
| 18 | +pfilesdir = $(prefix)/syspfiles | ||
| 19 | +pfiles_DATA = cterror.par |
src/cterror/cterror.cpp
0 → 100644
| 1 | +/*************************************************************************** | ||
| 2 | + * cterror - Parameter error calculation tool * | ||
| 3 | + * ----------------------------------------------------------------------- * | ||
| 4 | + * copyright (C) 2015 by Florent Forest * | ||
| 5 | + * ----------------------------------------------------------------------- * | ||
| 6 | + * * | ||
| 7 | + * This program is free software: you can redistribute it and/or modify * | ||
| 8 | + * it under the terms of the GNU General Public License as published by * | ||
| 9 | + * the Free Software Foundation, either version 3 of the License, or * | ||
| 10 | + * (at your option) any later version. * | ||
| 11 | + * * | ||
| 12 | + * This program is distributed in the hope that it will be useful, * | ||
| 13 | + * but WITHOUT ANY WARRANTY; without even the implied warranty of * | ||
| 14 | + * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the * | ||
| 15 | + * GNU General Public License for more details. * | ||
| 16 | + * * | ||
| 17 | + * You should have received a copy of the GNU General Public License * | ||
| 18 | + * along with this program. If not, see <http://www.gnu.org/licenses/>. * | ||
| 19 | + * * | ||
| 20 | + ***************************************************************************/ | ||
| 21 | +/** | ||
| 22 | + * @file cterror.cpp | ||
| 23 | + * @brief Parameter error calculation tool interface implementation | ||
| 24 | + * @author Florent Forest | ||
| 25 | + */ | ||
| 26 | + | ||
| 27 | +/* __ Includes ___________________________________________________________ */ | ||
| 28 | +#ifdef HAVE_CONFIG_H | ||
| 29 | +#include <config.h> | ||
| 30 | +#endif | ||
| 31 | +#include <stdio> | ||
| 32 | +#include "cterror.hpp" | ||
| 33 | +#include "GTools.hpp" | ||
| 34 | +#include "GOptimizer.hpp" | ||
| 35 | + | ||
| 36 | +/* __ Method name definitions ____________________________________________ */ | ||
| 37 | +#define G_GET_PARAMETERS "cterror::get_parameters()" | ||
| 38 | +#define G_GET_MODEL_PARAMETER "cterror::get_model_parameter()" | ||
| 39 | +#define G_EVALUATE "cterror::evaluate(double&)" | ||
| 40 | + | ||
| 41 | +/* __ Debug definitions __________________________________________________ */ | ||
| 42 | + | ||
| 43 | +/* __ Coding definitions _________________________________________________ */ | ||
| 44 | + | ||
| 45 | + | ||
| 46 | +/*========================================================================== | ||
| 47 | + = = | ||
| 48 | + = Constructors/destructors = | ||
| 49 | + = = | ||
| 50 | + ==========================================================================*/ | ||
| 51 | + | ||
| 52 | +/***********************************************************************//** | ||
| 53 | + * @brief Void constructor | ||
| 54 | + ***************************************************************************/ | ||
| 55 | +cterror::cterror(void) : ctool(CTULIMIT_NAME, CTULIMIT_VERSION) | ||
| 56 | +{ | ||
| 57 | + // Initialise members | ||
| 58 | + init_members(); | ||
| 59 | + | ||
| 60 | + // Return | ||
| 61 | + return; | ||
| 62 | +} | ||
| 63 | + | ||
| 64 | + | ||
| 65 | +/***********************************************************************//** | ||
| 66 | + * @brief Observations constructor | ||
| 67 | + * | ||
| 68 | + * param[in] obs Observation container. | ||
| 69 | + * | ||
| 70 | + * This method creates an instance of the class by copying an existing | ||
| 71 | + * observations container. | ||
| 72 | + ***************************************************************************/ | ||
| 73 | +cterror::cterror(const GObservations& obs) : | ||
| 74 | + ctool(CTULIMIT_NAME, CTULIMIT_VERSION) | ||
| 75 | +{ | ||
| 76 | + // Initialise members | ||
| 77 | + init_members(); | ||
| 78 | + | ||
| 79 | + // Set observations | ||
| 80 | + m_obs = obs; | ||
| 81 | + | ||
| 82 | + // Return | ||
| 83 | + return; | ||
| 84 | +} | ||
| 85 | + | ||
| 86 | + | ||
| 87 | + | ||
| 88 | +/***********************************************************************//** | ||
| 89 | + * @brief Command line constructor | ||
| 90 | + * | ||
| 91 | + * @param[in] argc Number of arguments in command line. | ||
| 92 | + * @param[in] argv Array of command line arguments. | ||
| 93 | + ***************************************************************************/ | ||
| 94 | +cterror::cterror(int argc, char *argv[]) : | ||
| 95 | + ctool(CTULIMIT_NAME, CTULIMIT_VERSION, argc, argv) | ||
| 96 | +{ | ||
| 97 | + // Initialise members | ||
| 98 | + init_members(); | ||
| 99 | + | ||
| 100 | + // Return | ||
| 101 | + return; | ||
| 102 | +} | ||
| 103 | + | ||
| 104 | + | ||
| 105 | +/***********************************************************************//** | ||
| 106 | + * @brief Copy constructor | ||
| 107 | + * | ||
| 108 | + * @param[in] app Application. | ||
| 109 | + ***************************************************************************/ | ||
| 110 | +cterror::cterror(const cterror& app) : ctool(app) | ||
| 111 | +{ | ||
| 112 | + // Initialise members | ||
| 113 | + init_members(); | ||
| 114 | + | ||
| 115 | + // Copy members | ||
| 116 | + copy_members(app); | ||
| 117 | + | ||
| 118 | + // Return | ||
| 119 | + return; | ||
| 120 | +} | ||
| 121 | + | ||
| 122 | + | ||
| 123 | +/***********************************************************************//** | ||
| 124 | + * @brief Destructor | ||
| 125 | + ***************************************************************************/ | ||
| 126 | +cterror::~cterror(void) | ||
| 127 | +{ | ||
| 128 | + // Free members | ||
| 129 | + free_members(); | ||
| 130 | + | ||
| 131 | + // Return | ||
| 132 | + return; | ||
| 133 | +} | ||
| 134 | + | ||
| 135 | + | ||
| 136 | +/*========================================================================== | ||
| 137 | + = = | ||
| 138 | + = Operators = | ||
| 139 | + = = | ||
| 140 | + ==========================================================================*/ | ||
| 141 | + | ||
| 142 | +/***********************************************************************//** | ||
| 143 | + * @brief Assignment operator | ||
| 144 | + * | ||
| 145 | + * @param[in] app Application. | ||
| 146 | + * @return Application. | ||
| 147 | + ***************************************************************************/ | ||
| 148 | +cterror& cterror::operator=(const cterror& app) | ||
| 149 | +{ | ||
| 150 | + // Execute only if object is not identical | ||
| 151 | + if (this != &app) { | ||
| 152 | + | ||
| 153 | + // Copy base class members | ||
| 154 | + this->ctool::operator=(app); | ||
| 155 | + | ||
| 156 | + // Free members | ||
| 157 | + free_members(); | ||
| 158 | + | ||
| 159 | + // Initialise members | ||
| 160 | + init_members(); | ||
| 161 | + | ||
| 162 | + // Copy members | ||
| 163 | + copy_members(app); | ||
| 164 | + | ||
| 165 | + } // endif: object was not identical | ||
| 166 | + | ||
| 167 | + // Return this object | ||
| 168 | + return *this; | ||
| 169 | +} | ||
| 170 | + | ||
| 171 | + | ||
| 172 | +/*========================================================================== | ||
| 173 | + = = | ||
| 174 | + = Public methods = | ||
| 175 | + = = | ||
| 176 | + ==========================================================================*/ | ||
| 177 | + | ||
| 178 | +/***********************************************************************//** | ||
| 179 | + * @brief Clear instance | ||
| 180 | + ***************************************************************************/ | ||
| 181 | +void cterror::clear(void) | ||
| 182 | +{ | ||
| 183 | + // Free members | ||
| 184 | + free_members(); | ||
| 185 | + this->ctool::free_members(); | ||
| 186 | + this->GApplication::free_members(); | ||
| 187 | + | ||
| 188 | + // Initialise members | ||
| 189 | + this->GApplication::init_members(); | ||
| 190 | + this->ctool::init_members(); | ||
| 191 | + init_members(); | ||
| 192 | + | ||
| 193 | + // Return | ||
| 194 | + return; | ||
| 195 | +} | ||
| 196 | + | ||
| 197 | + | ||
| 198 | +/***********************************************************************//** | ||
| 199 | + * @brief Computes | ||
| 200 | + * | ||
| 201 | + ***************************************************************************/ | ||
| 202 | +void cterror::run(void) | ||
| 203 | +{ | ||
| 204 | + // If we're in debug mode then all output is also dumped on the screen | ||
| 205 | + if (logDebug()) { | ||
| 206 | + log.cout(true); | ||
| 207 | + } | ||
| 208 | + | ||
| 209 | + // Get task parameters | ||
| 210 | + get_parameters(); | ||
| 211 | + | ||
| 212 | + // Write parameters into logger | ||
| 213 | + if (logTerse()) { | ||
| 214 | + log_parameters(); | ||
| 215 | + log << std::endl; | ||
| 216 | + } | ||
| 217 | + | ||
| 218 | + // Write observation(s) into logger | ||
| 219 | + if (logTerse()) { | ||
| 220 | + log << std::endl; | ||
| 221 | + if (m_obs.size() > 1) { | ||
| 222 | + log.header1("Observations"); | ||
| 223 | + } | ||
| 224 | + else { | ||
| 225 | + log.header1("Observation"); | ||
| 226 | + } | ||
| 227 | + log << m_obs << std::endl; | ||
| 228 | + } | ||
| 229 | + | ||
| 230 | + // Save original models | ||
| 231 | + GModels models_orig = m_obs.models(); | ||
| 232 | + | ||
| 233 | + // Save original log-likelihood. If the value is zero it has never been | ||
| 234 | + // computed hence we compute it now. | ||
| 235 | + m_best_logL = m_obs.logL(); | ||
| 236 | + if (m_best_logL == 0.0) { | ||
| 237 | + | ||
| 238 | + // Write header | ||
| 239 | + if (logTerse()) { | ||
| 240 | + log << std::endl; | ||
| 241 | + log.header1("Compute best-fit likelihood"); | ||
| 242 | + } | ||
| 243 | + | ||
| 244 | + // Reoptimize if likelihood was not given before | ||
| 245 | + GOptimizerLM* opt = new GOptimizerLM(); | ||
| 246 | + m_obs.optimize(*opt); | ||
| 247 | + m_obs.errors(*opt); | ||
| 248 | + m_best_logL = m_obs.logL(); | ||
| 249 | + | ||
| 250 | + // Write optimised model into logger | ||
| 251 | + if (logTerse()) { | ||
| 252 | + log << m_obs.models() << std::endl; | ||
| 253 | + } | ||
| 254 | + | ||
| 255 | + } // endif: likelihood was zero | ||
| 256 | + | ||
| 257 | + // Compute parameter bracketing | ||
| 258 | + double value = m_model_par->value(); | ||
| 259 | + double error = m_model_par->error(); | ||
| 260 | + double parmin = value + m_sigma_min * error; | ||
| 261 | + double parmax = value + m_sigma_max * error; | ||
| 262 | + | ||
| 263 | + // Write header | ||
| 264 | + if (logTerse()) { | ||
| 265 | + log << std::endl; | ||
| 266 | + log.header1("Compute upper limit"); | ||
| 267 | + log << gammalib::parformat("Model name"); | ||
| 268 | + log << m_skymodel->name() << std::endl; | ||
| 269 | + log << gammalib::parformat("Parameter name"); | ||
| 270 | + log << m_model_par->name() << std::endl; | ||
| 271 | + log << gammalib::parformat("Confidence level"); | ||
| 272 | + log << m_confidence*100.0 << "%" << std::endl; | ||
| 273 | + log << gammalib::parformat("Log-likelihood difference"); | ||
| 274 | + log << m_dlogL << std::endl; | ||
| 275 | + log << gammalib::parformat("Initial parameter range"); | ||
| 276 | + log << "["; | ||
| 277 | + log << parmin; | ||
| 278 | + log << ", "; | ||
| 279 | + log << parmax; | ||
| 280 | + log << "]" << std::endl; | ||
| 281 | + } | ||
| 282 | + | ||
| 283 | + // Compute lower boundary | ||
| 284 | + error_bisection(parmin, value); | ||
| 285 | + double value_lo = m_model_par->value(); | ||
| 286 | + | ||
| 287 | + if (logTerse()) { | ||
| 288 | + log << gammalib::parformat("Final lower parameter value"); | ||
| 289 | + log << value_lo << std::endl; | ||
| 290 | + } | ||
| 291 | + | ||
| 292 | + // Reset parameter to previous value | ||
| 293 | + m_model_par->value(value); | ||
| 294 | + | ||
| 295 | + // Compute upper boundary | ||
| 296 | + error_bisection(value, parmax); | ||
| 297 | + double value_hi = m_model_par->value(); | ||
| 298 | + | ||
| 299 | + if (logTerse()) { | ||
| 300 | + log << gammalib::parformat("Final upper parameter"); | ||
| 301 | + log << value_hi << std::endl; | ||
| 302 | + } | ||
| 303 | + | ||
| 304 | + // Compute error | ||
| 305 | + double error = value_hi - value_lo; | ||
| 306 | + | ||
| 307 | + // Write results to logfile | ||
| 308 | + if (logTerse()) { | ||
| 309 | + log << std::endl; | ||
| 310 | + log.header1("Error results"); | ||
| 311 | + log << gammalib::parformat("Parameter error"); | ||
| 312 | + log << error; | ||
| 313 | + log << " " << m_model_par->unit(); | ||
| 314 | + log << std::endl; | ||
| 315 | + } | ||
| 316 | + | ||
| 317 | + // Recover original models | ||
| 318 | + m_obs.models(models_orig); | ||
| 319 | + | ||
| 320 | + // Return | ||
| 321 | + return; | ||
| 322 | +} | ||
| 323 | + | ||
| 324 | +/*========================================================================== | ||
| 325 | + = = | ||
| 326 | + = Private methods = | ||
| 327 | + = = | ||
| 328 | + ==========================================================================*/ | ||
| 329 | + | ||
| 330 | +/***********************************************************************//** | ||
| 331 | + * @brief Initialise class members | ||
| 332 | + ***************************************************************************/ | ||
| 333 | +void cterror::init_members(void) | ||
| 334 | +{ | ||
| 335 | + // Initialise user parameters | ||
| 336 | + m_srcname.clear(); | ||
| 337 | + m_confidence = 0.95; | ||
| 338 | + m_sigma_min = 0.0; | ||
| 339 | + m_sigma_max = 0.0; | ||
| 340 | + m_tol = 1.0e-6; | ||
| 341 | + m_max_iter = 50; | ||
| 342 | + | ||
| 343 | + // Initialise protected members | ||
| 344 | + m_obs.clear(); | ||
| 345 | + m_dlogL = 0.0; | ||
| 346 | + m_skymodel = NULL; | ||
| 347 | + m_model_par = NULL; | ||
| 348 | + m_best_logL = 0.0; | ||
| 349 | + | ||
| 350 | + // Return | ||
| 351 | + return; | ||
| 352 | +} | ||
| 353 | + | ||
| 354 | + | ||
| 355 | +/***********************************************************************//** | ||
| 356 | + * @brief Copy class members | ||
| 357 | + * | ||
| 358 | + * @param[in] app Application. | ||
| 359 | + ***************************************************************************/ | ||
| 360 | +void cterror::copy_members(const cterror& app) | ||
| 361 | +{ | ||
| 362 | + // Copy user parameters | ||
| 363 | + m_srcname = app.m_srcname; | ||
| 364 | + m_confidence = app.m_confidence; | ||
| 365 | + m_sigma_min = app.m_sigma_min; | ||
| 366 | + m_sigma_max = app.m_sigma_max; | ||
| 367 | + m_tol = app.m_tol; | ||
| 368 | + m_max_iter = app.m_max_iter; | ||
| 369 | + | ||
| 370 | + | ||
| 371 | + // Copy protected members | ||
| 372 | + m_obs = app.m_obs; | ||
| 373 | + m_dlogL = app.m_dlogL; | ||
| 374 | + m_best_logL = app.m_best_logL; | ||
| 375 | + | ||
| 376 | + | ||
| 377 | + // Extract model parameter | ||
| 378 | + get_model_parameter(); | ||
| 379 | + | ||
| 380 | + // Return | ||
| 381 | + return; | ||
| 382 | +} | ||
| 383 | + | ||
| 384 | + | ||
| 385 | +/***********************************************************************//** | ||
| 386 | + * @brief Delete class members | ||
| 387 | + ***************************************************************************/ | ||
| 388 | +void cterror::free_members(void) | ||
| 389 | +{ | ||
| 390 | + // Return | ||
| 391 | + return; | ||
| 392 | +} | ||
| 393 | + | ||
| 394 | + | ||
| 395 | +/***********************************************************************//** | ||
| 396 | + * @brief Get application parameters | ||
| 397 | + * | ||
| 398 | + * @exception GException::invalid_value | ||
| 399 | + * Test source not found. | ||
| 400 | + * | ||
| 401 | + * Get all task parameters from parameter file or (if required) by querying | ||
| 402 | + * the user. | ||
| 403 | + ***************************************************************************/ | ||
| 404 | +void cterror::get_parameters(void) | ||
| 405 | +{ | ||
| 406 | + // If there are no observations in container then load them via user | ||
| 407 | + // parameters | ||
| 408 | + if (m_obs.size() == 0) { | ||
| 409 | + | ||
| 410 | + // Throw exception if no input observation file is given | ||
| 411 | + require_inobs(G_GET_PARAMETERS); | ||
| 412 | + | ||
| 413 | + // Build observation container | ||
| 414 | + m_obs = get_observations(); | ||
| 415 | + | ||
| 416 | + } // endif: there was no observation in the container | ||
| 417 | + | ||
| 418 | + | ||
| 419 | + // If there is are no models associated with the observations then | ||
| 420 | + // load now the model definition | ||
| 421 | + if (m_obs.models().size() == 0) { | ||
| 422 | + | ||
| 423 | + // Get models XML filename | ||
| 424 | + std::string filename = (*this)["inmodel"].filename(); | ||
| 425 | + | ||
| 426 | + // Setup models for optimizing. | ||
| 427 | + m_obs.models(GModels(filename)); | ||
| 428 | + | ||
| 429 | + } // endif: no models were associated with observations | ||
| 430 | + | ||
| 431 | + // Get name of test source and check container for this name | ||
| 432 | + m_srcname = (*this)["srcname"].string(); | ||
| 433 | + if (!m_obs.models().contains(m_srcname)) { | ||
| 434 | + std::string msg = "Source \""+m_srcname+"\" not found in model " | ||
| 435 | + "container. Please add a source with that name " | ||
| 436 | + "or check for possible typos."; | ||
| 437 | + throw GException::invalid_value(G_GET_PARAMETERS, msg); | ||
| 438 | + } | ||
| 439 | + | ||
| 440 | + // Get relevant model and parameter for upper limit computation | ||
| 441 | + get_model_parameter(); | ||
| 442 | + | ||
| 443 | + // Get confidence level and transform into log-likelihood difference | ||
| 444 | + m_confidence = (*this)["confidence"].real(); | ||
| 445 | + double sigma = gammalib::erfinv(m_confidence) * gammalib::sqrt_two; | ||
| 446 | + m_dlogL = (sigma*sigma) / 2.0; | ||
| 447 | + | ||
| 448 | + // Read starting boundaries for bisection | ||
| 449 | + m_sigma_min = (*this)["sigma_min"].real(); | ||
| 450 | + m_sigma_max = (*this)["sigma_max"].real(); | ||
| 451 | + | ||
| 452 | + // Read energy values | ||
| 453 | + m_eref = (*this)["eref"].real(); | ||
| 454 | + m_emin = (*this)["emin"].real(); | ||
| 455 | + m_emax = (*this)["emax"].real(); | ||
| 456 | + | ||
| 457 | + // Read precision | ||
| 458 | + m_tol = (*this)["tol"].real(); | ||
| 459 | + m_max_iter = (*this)["max_iter"].integer(); | ||
| 460 | + | ||
| 461 | + // Return | ||
| 462 | + return; | ||
| 463 | +} | ||
| 464 | + | ||
| 465 | + | ||
| 466 | +/***********************************************************************//** | ||
| 467 | + * @brief Get application parameters | ||
| 468 | + * | ||
| 469 | + * @exception GException::invalid_value | ||
| 470 | + * Did not find a valid model parameter | ||
| 471 | + * | ||
| 472 | + * Extracts a pointer to the sky model (m_skymodel) and a pointer to the | ||
| 473 | + * relevant model parameter (m_model_par) from the model container. | ||
| 474 | + ***************************************************************************/ | ||
| 475 | +void cterror::get_model_parameter(void) | ||
| 476 | +{ | ||
| 477 | + // Continue only if source model exists | ||
| 478 | + if (m_obs.models().contains(m_srcname)) { | ||
| 479 | + | ||
| 480 | + // Get relevant model and parameter for upper limit computation. | ||
| 481 | + GModels& models = const_cast<GModels&>(m_obs.models()); | ||
| 482 | + m_skymodel = dynamic_cast<GModelSky*>(models[m_srcname]); | ||
| 483 | + if (m_skymodel == NULL) { | ||
| 484 | + std::string msg = "Source \""+m_srcname+"\" is not a sky model. " | ||
| 485 | + "Please specify the name of a sky model for " | ||
| 486 | + "upper limit computation."; | ||
| 487 | + throw GException::invalid_value(G_GET_MODEL_PARAMETER, msg); | ||
| 488 | + } | ||
| 489 | + if (m_skymodel->spectral()->type() == "NodeFunction") { | ||
| 490 | + std::string msg = "\"NodeFunction\" cannot be used as spectral " | ||
| 491 | + "model for an upper limit computation. " | ||
| 492 | + "Please specify another spectral model."; | ||
| 493 | + throw GException::invalid_value(G_GET_MODEL_PARAMETER, msg); | ||
| 494 | + } | ||
| 495 | + if (m_skymodel->spectral()->has_par("Normalization")) { | ||
| 496 | + m_model_par = &(m_skymodel->spectral()->operator[]("Normalization")); | ||
| 497 | + } | ||
| 498 | + else if (m_skymodel->spectral()->has_par("Prefactor")) { | ||
| 499 | + m_model_par = &(m_skymodel->spectral()->operator[]("Prefactor")); | ||
| 500 | + } | ||
| 501 | + else if (m_skymodel->spectral()->has_par("Integral")) { | ||
| 502 | + m_model_par = &(m_skymodel->spectral()->operator[]("Integral")); | ||
| 503 | + } | ||
| 504 | + else if (m_skymodel->spectral()->has_par("Value")) { | ||
| 505 | + m_model_par = &(m_skymodel->spectral()->operator[]("Value")); | ||
| 506 | + } | ||
| 507 | + else { | ||
| 508 | + std::string msg = "Require spectral parameter \"Normalization\", " | ||
| 509 | + "\"Prefactor\", \"Integral\" or \"Value\" for " | ||
| 510 | + "upper limit computation. The specified source " | ||
| 511 | + "\""+m_srcname+"\" does not have such a parameter."; | ||
| 512 | + throw GException::invalid_value(G_GET_MODEL_PARAMETER, msg); | ||
| 513 | + } | ||
| 514 | + | ||
| 515 | + } // endif: source model existed | ||
| 516 | + | ||
| 517 | + // Return | ||
| 518 | + return; | ||
| 519 | +} | ||
| 520 | + | ||
| 521 | + | ||
| 522 | +/***********************************************************************//** | ||
| 523 | + * @brief Performs error computation by using a bisection method | ||
| 524 | + * | ||
| 525 | + * @param[in] min Minimum parameter value | ||
| 526 | + * @param[in] max Maximum parameter value | ||
| 527 | + * | ||
| 528 | + * This method calculates the error using a bisection method. | ||
| 529 | + ***************************************************************************/ | ||
| 530 | +void cterror::error_bisection(const double& min, const double& max) | ||
| 531 | +{ | ||
| 532 | + // Copy values to working values | ||
| 533 | + double wrk_min = min; | ||
| 534 | + double wrk_max = max; | ||
| 535 | + | ||
| 536 | + // Initialise iteration counter | ||
| 537 | + int iter = 0; | ||
| 538 | + | ||
| 539 | + // Loop until breaking condition is reached | ||
| 540 | + while (true) { | ||
| 541 | + | ||
| 542 | + // Log information | ||
| 543 | + if (logExplicit()) { | ||
| 544 | + log << gammalib::parformat("Iteration "+gammalib::str(iter)); | ||
| 545 | + log << "["; | ||
| 546 | + log << wrk_min; | ||
| 547 | + log << ", "; | ||
| 548 | + log << wrk_max; | ||
| 549 | + log << "]" << std::endl; | ||
| 550 | + } | ||
| 551 | + | ||
| 552 | + // Throw exception if maximum iterations are reached | ||
| 553 | + if (iter > m_max_iter) { | ||
| 554 | + std::string msg = "The maximum number of "+gammalib::str(m_max_iter)+ | ||
| 555 | + " has been reached. You may consider to increase" | ||
| 556 | + " the \"max_iter\" parameter and re-run cterror."; | ||
| 557 | + throw GException::invalid_value(G_UL_BISECTION, msg); | ||
| 558 | + } | ||
| 559 | + | ||
| 560 | + // Compute center of boundary | ||
| 561 | + double mid = (wrk_min + wrk_max) / 2.0; | ||
| 562 | + | ||
| 563 | + // Calculate function value | ||
| 564 | + double eval_mid = evaluate(mid); | ||
| 565 | + | ||
| 566 | + // Check for convergence inside tolerance | ||
| 567 | + if (std::abs(eval_mid) < m_tol) { | ||
| 568 | + break; | ||
| 569 | + } | ||
| 570 | + | ||
| 571 | + // Change boundaries for further iteration | ||
| 572 | + if (eval_mid > 0.0) { | ||
| 573 | + wrk_max = mid; | ||
| 574 | + } | ||
| 575 | + else if (eval_mid < 0.0) { | ||
| 576 | + wrk_min = mid; | ||
| 577 | + } | ||
| 578 | + | ||
| 579 | + // Increment counter | ||
| 580 | + iter++; | ||
| 581 | + | ||
| 582 | + } // endwhile | ||
| 583 | + | ||
| 584 | + // Return | ||
| 585 | + return; | ||
| 586 | + | ||
| 587 | +} | ||
| 588 | + | ||
| 589 | + | ||
| 590 | +/***********************************************************************//** | ||
| 591 | + * @brief Evaluates the log-likelihood | ||
| 592 | + * | ||
| 593 | + * @param[in] value Parameter factor value | ||
| 594 | + * @return Log-likelihood value | ||
| 595 | + * | ||
| 596 | + * This method evaluates the log-likelihood as a function of the parameter | ||
| 597 | + * of interest. | ||
| 598 | + ***************************************************************************/ | ||
| 599 | +double cterror::evaluate(const double& value) | ||
| 600 | +{ | ||
| 601 | + // Initialise log-likelihood value | ||
| 602 | + double logL = 0.0; | ||
| 603 | + | ||
| 604 | + // Check if given parameter is within boundaries | ||
| 605 | + if (value > m_model_par->min() && value < m_model_par->max()) { | ||
| 606 | + | ||
| 607 | + // Change parameter factor | ||
| 608 | + m_model_par->value(value); | ||
| 609 | + | ||
| 610 | + // Evaluate likelihood for new model container | ||
| 611 | + m_obs.eval(); | ||
| 612 | + | ||
| 613 | + // Retrieve likelihood | ||
| 614 | + logL = m_obs.logL(); | ||
| 615 | + | ||
| 616 | + } // endif: value was inside allowed range | ||
| 617 | + | ||
| 618 | + // ... otherwise signal that the parameter went outside the boundaries | ||
| 619 | + else { | ||
| 620 | + std::string msg = "Value of parameter \""+m_model_par->name()+"\" " | ||
| 621 | + "outside of validity range requested. To omit " | ||
| 622 | + "this error please enlarge the parameter range " | ||
| 623 | + "in the model XML file."; | ||
| 624 | + throw GException::invalid_value(G_EVALUATE, msg); | ||
| 625 | + } | ||
| 626 | + | ||
| 627 | + // Compute function value | ||
| 628 | + double logL_difference = logL - m_best_logL - m_dlogL; | ||
| 629 | + | ||
| 630 | + // Return | ||
| 631 | + return logL_difference; | ||
| 632 | +} |
src/cterror/cterror.hpp
0 → 100644
| 1 | +/*************************************************************************** | ||
| 2 | + * cterror - Parameter error calculation tool * | ||
| 3 | + * ----------------------------------------------------------------------- * | ||
| 4 | + * copyright (C) 2015 by Florent Forest * | ||
| 5 | + * ----------------------------------------------------------------------- * | ||
| 6 | + * * | ||
| 7 | + * This program is free software: you can redistribute it and/or modify * | ||
| 8 | + * it under the terms of the GNU General Public License as published by * | ||
| 9 | + * the Free Software Foundation, either version 3 of the License, or * | ||
| 10 | + * (at your option) any later version. * | ||
| 11 | + * * | ||
| 12 | + * This program is distributed in the hope that it will be useful, * | ||
| 13 | + * but WITHOUT ANY WARRANTY; without even the implied warranty of * | ||
| 14 | + * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the * | ||
| 15 | + * GNU General Public License for more details. * | ||
| 16 | + * * | ||
| 17 | + * You should have received a copy of the GNU General Public License * | ||
| 18 | + * along with this program. If not, see <http://www.gnu.org/licenses/>. * | ||
| 19 | + * * | ||
| 20 | + ***************************************************************************/ | ||
| 21 | +/** | ||
| 22 | + * @file cterror.hpp | ||
| 23 | + * @brief Parameter error calculation tool interface definition | ||
| 24 | + * @author Florent Forest | ||
| 25 | + */ | ||
| 26 | + | ||
| 27 | +#ifndef CTERROR_HPP | ||
| 28 | +#define CTERROR_HPP | ||
| 29 | + | ||
| 30 | +/* __ Includes ___________________________________________________________ */ | ||
| 31 | +#include "GammaLib.hpp" | ||
| 32 | +#include "GCTALib.hpp" | ||
| 33 | +#include "ctool.hpp" | ||
| 34 | + | ||
| 35 | +/* __Definitions _________________________________________________________ */ | ||
| 36 | +#define CTERROR_NAME "cterror" | ||
| 37 | +#define CTERROR_VERSION "01-00-00" | ||
| 38 | + | ||
| 39 | + | ||
| 40 | +/***********************************************************************//** | ||
| 41 | + * @class cterror | ||
| 42 | + * | ||
| 43 | + * @brief Parameter error calculation tool | ||
| 44 | + * | ||
| 45 | + ***************************************************************************/ | ||
| 46 | +class cterror : public ctool { | ||
| 47 | + | ||
| 48 | +public: | ||
| 49 | + // Constructors and destructors | ||
| 50 | + cterror(void); | ||
| 51 | + explicit cterror(const GObservations& obs); | ||
| 52 | + cterror(int argc, char *argv[]); | ||
| 53 | + cterror(const cterror& app); | ||
| 54 | + virtual ~cterror(void); | ||
| 55 | + | ||
| 56 | + // Operators | ||
| 57 | + cterror& operator=(const cterror& app); | ||
| 58 | + | ||
| 59 | + // Methods | ||
| 60 | + void clear(void); | ||
| 61 | + void run(void); | ||
| 62 | + void save(void); | ||
| 63 | + const GObservations& obs(void) const; | ||
| 64 | + | ||
| 65 | + | ||
| 66 | +protected: | ||
| 67 | + // Protected methods | ||
| 68 | + void init_members(void); | ||
| 69 | + void copy_members(const cterror& app); | ||
| 70 | + void free_members(void); | ||
| 71 | + void get_parameters(void); | ||
| 72 | + void get_model_parameter(void); | ||
| 73 | + double evaluate(const double& value); | ||
| 74 | + void error_bisection(const double& min, const double& max); | ||
| 75 | + | ||
| 76 | + // User parameters | ||
| 77 | + std::string m_srcname; //!< Name of source which is moved around | ||
| 78 | + double m_confidence; //!< Confidence level | ||
| 79 | + double m_sigma_min; //!< Starting value minimum (multiple fit errors above fit values) | ||
| 80 | + double m_sigma_max; //!< Starting value maximum (multiple fit errors above fit values) | ||
| 81 | + double m_tol; //!< Tolerance for limit determination | ||
| 82 | + int m_max_iter; //!< Maximum number of iterations | ||
| 83 | + | ||
| 84 | + // Protected members | ||
| 85 | + GObservations m_obs; //!< Observation container | ||
| 86 | + double m_dlogL; //!< Likelihood difference for upper limit computation | ||
| 87 | + GModelSky* m_skymodel; //!< Pointer to sky model | ||
| 88 | + GModelPar* m_model_par; //!< Pointer to model parameter | ||
| 89 | + double m_best_logL; //!< Best fit log likelihood of given model | ||
| 90 | + | ||
| 91 | +}; | ||
| 92 | + | ||
| 93 | + | ||
| 94 | +/***********************************************************************//** | ||
| 95 | + * @brief Return observation container | ||
| 96 | + * | ||
| 97 | + * @return Reference to observation container | ||
| 98 | + ***************************************************************************/ | ||
| 99 | +inline | ||
| 100 | +const GObservations& cterror::obs(void) const | ||
| 101 | +{ | ||
| 102 | + return m_obs; | ||
| 103 | +} | ||
| 104 | + | ||
| 105 | + | ||
| 106 | +#endif /* CTERROR_HPP */ |
src/cterror/cterror.par
0 → 100644
| 1 | +############################################################################# | ||
| 2 | +# # | ||
| 3 | +# CTA Science Analysis Tools # | ||
| 4 | +# # | ||
| 5 | +# ------------------------------------------------------------------------- # | ||
| 6 | +# # | ||
| 7 | +# File: cterror.par # | ||
| 8 | +# Version: 01-00-00 # | ||
| 9 | +# # | ||
| 10 | +# Author: Florent Forest # | ||
| 11 | +# florent.forest@isae.fr # | ||
| 12 | +# ISAE-Supaero # | ||
| 13 | +# # | ||
| 14 | +# Purpose: Parameter file for the cterror application # | ||
| 15 | +# # | ||
| 16 | +############################################################################# | ||
| 17 | +# | ||
| 18 | +# General parameters | ||
| 19 | +#=================== | ||
| 20 | +inobs, f, a, "events.fits",,, "Input event list, counts cube or observation definition file" | ||
| 21 | +inmodel, f, a, "$CTOOLS/share/models/crab.xml",,, "Source model" | ||
| 22 | +srcname, s, a, "Crab",,, "Source of interest" | ||
| 23 | +expcube, f, a, "NONE",,, "Exposure cube file (only needed for stacked analysis)" | ||
| 24 | +psfcube, f, a, "NONE",,, "PSF cube file (only needed for stacked analysis)" | ||
| 25 | +bkgcube, f, a, "NONE",,, "Background cube file (only needed for stacked analysis)" | ||
| 26 | +caldb, s, a, "dummy",,, "Calibration database" | ||
| 27 | +irf, s, a, "cta_dummy_irf",,, "Instrument response function" | ||
| 28 | + | ||
| 29 | +# | ||
| 30 | +# Upper limit calculation parameters | ||
| 31 | +#=================================== | ||
| 32 | +confidence, r, h, 0.68,0.0,1.0, "Confidence level" | ||
| 33 | +sigma_min, r, h, 0.0,,, "Start minimum value (multiple sigmas above fit value)" | ||
| 34 | +sigma_max, r, h, 10.0,,, "Start maximum value (multiple sigmas above fit value)" | ||
| 35 | +eref, r, h, 1.0,,, "Reference energy for differential limits (TeV)" | ||
| 36 | +emin, r, h, 1.0,,, "Minimum energy for flux limits (TeV)" | ||
| 37 | +emax, r, h, 100.0,,, "Maximum energy for flux limits (TeV)" | ||
| 38 | +tol, r, h, 1e-5,,, "Computation tolerance" | ||
| 39 | +max_iter, i, h, 50,,, "Maximum number of iterations" | ||
| 40 | + | ||
| 41 | +# | ||
| 42 | +# Standard parameters | ||
| 43 | +#==================== | ||
| 44 | +chatter, i, h, 2,0,4, "Chattiness of output" | ||
| 45 | +clobber, b, h, yes,,, "Overwrite existing output files with new output files?" | ||
| 46 | +debug, b, h, no,,, "Debugging mode activated" | ||
| 47 | +mode, s, h, "ql",,, "Mode of automatic parameters" | ||
| 48 | +logfile, f, h, "cterror.log",,, "Log filename" |
src/cterror/main.cpp
0 → 100644
| 1 | +/*************************************************************************** | ||
| 2 | + * cterror - Parameter error calculation tool * | ||
| 3 | + * ----------------------------------------------------------------------- * | ||
| 4 | + * copyright (C) 2015 by Florent Forest * | ||
| 5 | + * ----------------------------------------------------------------------- * | ||
| 6 | + * * | ||
| 7 | + * This program is free software: you can redistribute it and/or modify * | ||
| 8 | + * it under the terms of the GNU General Public License as published by * | ||
| 9 | + * the Free Software Foundation, either version 3 of the License, or * | ||
| 10 | + * (at your option) any later version. * | ||
| 11 | + * * | ||
| 12 | + * This program is distributed in the hope that it will be useful, * | ||
| 13 | + * but WITHOUT ANY WARRANTY; without even the implied warranty of * | ||
| 14 | + * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the * | ||
| 15 | + * GNU General Public License for more details. * | ||
| 16 | + * * | ||
| 17 | + * You should have received a copy of the GNU General Public License * | ||
| 18 | + * along with this program. If not, see <http://www.gnu.org/licenses/>. * | ||
| 19 | + * * | ||
| 20 | + ***************************************************************************/ | ||
| 21 | +/** | ||
| 22 | + * @file main.cpp | ||
| 23 | + * @brief CTA parameter error calculation | ||
| 24 | + * @author Florent Forest | ||
| 25 | + */ | ||
| 26 | + | ||
| 27 | +/* __ Includes ___________________________________________________________ */ | ||
| 28 | +#ifdef HAVE_CONFIG_H | ||
| 29 | +#include <config.h> | ||
| 30 | +#endif | ||
| 31 | +#include <cstdio> | ||
| 32 | +#include "cterror.hpp" | ||
| 33 | + | ||
| 34 | + | ||
| 35 | +/***********************************************************************//** | ||
| 36 | + * @brief Main entry point of application | ||
| 37 | + * | ||
| 38 | + * @param[in] argc Number of command line arguments. | ||
| 39 | + * @param[in] argv Command line arguments. | ||
| 40 | + * | ||
| 41 | + * This is the main entry point of the cterror application. It allocates a | ||
| 42 | + * cterror object and runs the application. | ||
| 43 | + ***************************************************************************/ | ||
| 44 | +int main (int argc, char *argv[]) | ||
| 45 | +{ | ||
| 46 | + // Initialise return code | ||
| 47 | + int rc = 1; | ||
| 48 | + | ||
| 49 | + // Create instance of application | ||
| 50 | + cterror application(argc, argv); | ||
| 51 | + | ||
| 52 | + // Run application | ||
| 53 | + try { | ||
| 54 | + // Execute application | ||
| 55 | + application.execute(); | ||
| 56 | + | ||
| 57 | + // Signal success | ||
| 58 | + rc = 0; | ||
| 59 | + } | ||
| 60 | + catch (std::exception &e) { | ||
| 61 | + | ||
| 62 | + // Extract error message | ||
| 63 | + std::string message = e.what(); | ||
| 64 | + std::string signal = "*** ERROR encounterted in the execution of" | ||
| 65 | + " cterror. Run aborted ..."; | ||
| 66 | + | ||
| 67 | + // Write error in logger | ||
| 68 | + application.log << signal << std::endl; | ||
| 69 | + application.log << message << std::endl; | ||
| 70 | + | ||
| 71 | + // Write error on standard output | ||
| 72 | + std::cout << signal << std::endl; | ||
| 73 | + std::cout << message << std::endl; | ||
| 74 | + | ||
| 75 | + } // endcatch: catched any application error | ||
| 76 | + | ||
| 77 | + // Return | ||
| 78 | + return rc; | ||
| 79 | +} |