Commit 4dc88fd65f7ae20f164316b4ad27f5f9296b9588
1 parent
996c4df2
Merge in cterror application.
Showing
9 changed files
with
142 additions
and
117 deletions
pyext/cterror.i
| 1 | 1 | /*************************************************************************** |
| 2 | - * cterror - Parameter error calculation tool * | |
| 2 | + * cterror - Parameter error calculation tool * | |
| 3 | 3 | * ----------------------------------------------------------------------- * |
| 4 | 4 | * copyright (C) 2015 by Florent Forest * |
| 5 | 5 | * ----------------------------------------------------------------------- * |
| ... | ... | @@ -38,7 +38,7 @@ |
| 38 | 38 | ***************************************************************************/ |
| 39 | 39 | class cterror : public ctool { |
| 40 | 40 | |
| 41 | - public: | |
| 41 | +public: | |
| 42 | 42 | // Constructors and destructors |
| 43 | 43 | cterror(void); |
| 44 | 44 | explicit cterror(const GObservations& obs); |
| ... | ... | @@ -52,9 +52,9 @@ class cterror : public ctool { |
| 52 | 52 | void save(void); |
| 53 | 53 | const GObservations& obs(void) const; |
| 54 | 54 | const GOptimizer* opt(void) const; |
| 55 | - | |
| 56 | 55 | }; |
| 57 | 56 | |
| 57 | + | |
| 58 | 58 | /***********************************************************************//** |
| 59 | 59 | * @brief Parameter error calculation tool Python extensions |
| 60 | 60 | ***************************************************************************/ | ... | ... |
pyext/ctools/obsutils.py
| ... | ... | @@ -140,9 +140,10 @@ def fit(obs, log=False, debug=False, chatter=2, edisp=False): |
| 140 | 140 | Parameters: |
| 141 | 141 | obs - Observation container |
| 142 | 142 | Keywords: |
| 143 | - log - Create log file(s) | |
| 144 | - debug - Create screen dump | |
| 145 | - edisp - Apply energy dispersion? | |
| 143 | + log - Create log file(s) | |
| 144 | + debug - Create screen dump | |
| 145 | + chatter - Chatter level | |
| 146 | + edisp - Apply energy dispersion? | |
| 146 | 147 | """ |
| 147 | 148 | # Allocate ctlike application |
| 148 | 149 | like = ctools.ctlike(obs) |
| ... | ... | @@ -167,19 +168,21 @@ def fit(obs, log=False, debug=False, chatter=2, edisp=False): |
| 167 | 168 | # Return observations |
| 168 | 169 | return like |
| 169 | 170 | |
| 170 | -# ================ # | |
| 171 | -# Run cterror # | |
| 172 | -# ================ # | |
| 173 | -def cterror(obs, log=False, debug=False): | |
| 171 | + | |
| 172 | +# ============================================================== # | |
| 173 | +# Fit observations and determine errors using likelihood profile # | |
| 174 | +# ============================================================== # | |
| 175 | +def cterror(obs, log=False, debug=False, chatter=2): | |
| 174 | 176 | """ |
| 175 | 177 | Perform maximum likelihood fitting of observations in the container. |
| 176 | 178 | |
| 177 | 179 | Parameters: |
| 178 | 180 | obs - Observation container |
| 179 | 181 | Keywords: |
| 180 | - log - Create log file(s) | |
| 181 | - debug - Create screen dump | |
| 182 | - edisp - Apply energy dispersion? | |
| 182 | + log - Create log file(s) | |
| 183 | + debug - Create screen dump | |
| 184 | + chatter - Chatter level | |
| 185 | + edisp - Apply energy dispersion? | |
| 183 | 186 | """ |
| 184 | 187 | # Allocate cterror application |
| 185 | 188 | error = ctools.cterror(obs) |
| ... | ... | @@ -192,6 +195,9 @@ def cterror(obs, log=False, debug=False): |
| 192 | 195 | if debug: |
| 193 | 196 | error["debug"].boolean(True) |
| 194 | 197 | |
| 198 | + # Set chatter level | |
| 199 | + error["chatter"].chatter = chatter | |
| 200 | + | |
| 195 | 201 | # Run cterror application. |
| 196 | 202 | error.run() |
| 197 | 203 | ... | ... |
scripts/cspull.py
| ... | ... | @@ -47,6 +47,8 @@ class cspull(ctools.cscript): |
| 47 | 47 | self.obs = None |
| 48 | 48 | self.model = None |
| 49 | 49 | self.m_inmodel = None |
| 50 | + self.m_edisp = False | |
| 51 | + self.m_profile = False | |
| 50 | 52 | |
| 51 | 53 | # Make sure that parfile exists |
| 52 | 54 | file = self.parfile() |
| ... | ... | @@ -103,6 +105,7 @@ class cspull(ctools.cscript): |
| 103 | 105 | pars.append(gammalib.GApplicationPar("irf","s","a","South_50h","","","Instrument response function")) |
| 104 | 106 | pars.append(gammalib.GApplicationPar("deadc","r","h","0.95","0","1","Deadtime correction factor")) |
| 105 | 107 | pars.append(gammalib.GApplicationPar("edisp","b","h","no","","","Apply energy dispersion?")) |
| 108 | + pars.append(gammalib.GApplicationPar("profile","b","h","no","","","Use likelihood profile method for errors?")) | |
| 106 | 109 | pars.append(gammalib.GApplicationPar("ntrials","i","a","10","","","Number of trials")) |
| 107 | 110 | pars.append(gammalib.GApplicationPar("ra","r","a","83.6331","0","360","RA of pointing (deg)")) |
| 108 | 111 | pars.append(gammalib.GApplicationPar("dec","r","a","22.0145","-90","90","Dec of pointing (deg)")) |
| ... | ... | @@ -160,6 +163,7 @@ class cspull(ctools.cscript): |
| 160 | 163 | self.m_ntrials = self["ntrials"].integer() |
| 161 | 164 | self.m_edisp = self["edisp"].boolean() |
| 162 | 165 | self.m_offset = self["offset"].real() |
| 166 | + self.m_profile = self["profile"].boolean() | |
| 163 | 167 | |
| 164 | 168 | # Set some fixed parameters |
| 165 | 169 | self.m_log = False # Logging in client tools |
| ... | ... | @@ -297,9 +301,12 @@ class cspull(ctools.cscript): |
| 297 | 301 | self.log("\n") |
| 298 | 302 | |
| 299 | 303 | # Fit model |
| 300 | - #like = obsutils.fit(obs, edisp=self.m_edisp, \ | |
| 301 | - # log=self.m_log, debug=self.m_debug) | |
| 302 | - like = obsutils.cterror(obs, log=self.m_log, debug=self.m_debug) | |
| 304 | + if self.m_profile: | |
| 305 | + like = obsutils.cterror(obs, \ | |
| 306 | + log=self.m_log, debug=self.m_debug) | |
| 307 | + else: | |
| 308 | + like = obsutils.fit(obs, edisp=self.m_edisp, \ | |
| 309 | + log=self.m_log, debug=self.m_debug) | |
| 303 | 310 | |
| 304 | 311 | # Store results |
| 305 | 312 | logL = like.opt().value() | ... | ... |
src/cterror/cterror.cpp
| 1 | 1 | /*************************************************************************** |
| 2 | - * cterror - Parameter error calculation tool * | |
| 2 | + * cterror - Parameter error calculation tool * | |
| 3 | 3 | * ----------------------------------------------------------------------- * |
| 4 | 4 | * copyright (C) 2015 by Florent Forest * |
| 5 | 5 | * ----------------------------------------------------------------------- * |
| ... | ... | @@ -72,7 +72,7 @@ cterror::cterror(void) : ctool(CTERROR_NAME, CTERROR_VERSION) |
| 72 | 72 | * observations container. |
| 73 | 73 | ***************************************************************************/ |
| 74 | 74 | cterror::cterror(const GObservations& obs) : |
| 75 | - ctool(CTERROR_NAME, CTERROR_VERSION) | |
| 75 | + ctool(CTERROR_NAME, CTERROR_VERSION) | |
| 76 | 76 | { |
| 77 | 77 | // Initialise members |
| 78 | 78 | init_members(); |
| ... | ... | @@ -93,7 +93,7 @@ cterror::cterror(const GObservations& obs) : |
| 93 | 93 | * @param[in] argv Array of command line arguments. |
| 94 | 94 | ***************************************************************************/ |
| 95 | 95 | cterror::cterror(int argc, char *argv[]) : |
| 96 | - ctool(CTERROR_NAME, CTERROR_VERSION, argc, argv) | |
| 96 | + ctool(CTERROR_NAME, CTERROR_VERSION, argc, argv) | |
| 97 | 97 | { |
| 98 | 98 | // Initialise members |
| 99 | 99 | init_members(); |
| ... | ... | @@ -198,7 +198,6 @@ void cterror::clear(void) |
| 198 | 198 | |
| 199 | 199 | /***********************************************************************//** |
| 200 | 200 | * @brief Computes |
| 201 | - * | |
| 202 | 201 | ***************************************************************************/ |
| 203 | 202 | void cterror::run(void) |
| 204 | 203 | { |
| ... | ... | @@ -228,11 +227,9 @@ void cterror::run(void) |
| 228 | 227 | log << m_obs << std::endl; |
| 229 | 228 | } |
| 230 | 229 | |
| 231 | - // Save original models | |
| 232 | - GModels models_orig = m_obs.models(); | |
| 233 | - | |
| 234 | 230 | // Allocate optimizer |
| 235 | 231 | m_opt = new GOptimizerLM(); |
| 232 | + | |
| 236 | 233 | // Optimize and save best log-likelihood |
| 237 | 234 | m_obs.optimize(*m_opt); |
| 238 | 235 | m_obs.errors(*m_opt); |
| ... | ... | @@ -249,35 +246,44 @@ void cterror::run(void) |
| 249 | 246 | log << m_obs.models() << std::endl; |
| 250 | 247 | } |
| 251 | 248 | |
| 252 | - | |
| 253 | 249 | // Continue only if source model exists |
| 254 | 250 | if (m_obs.models().contains(m_srcname)) { |
| 255 | 251 | |
| 256 | - GModels& models = const_cast<GModels&>(m_obs.models()); | |
| 257 | - m_skymodel = dynamic_cast<GModelSky*>(models[m_srcname]); | |
| 258 | - if (m_skymodel == NULL) { | |
| 252 | + // Save best fitting models | |
| 253 | + GModels models_best = m_obs.models(); | |
| 254 | + | |
| 255 | + // Get pointer on sky model | |
| 256 | + GModelSky* model = dynamic_cast<GModelSky*>(models_best[m_srcname]); | |
| 257 | + if (model == NULL) { | |
| 259 | 258 | std::string msg = "Source \""+m_srcname+"\" is not a sky model. " |
| 260 | - "Please specify the name of a sky model for " | |
| 261 | - "parameter error computation."; | |
| 259 | + "Please specify the name of a sky model for " | |
| 260 | + "parameter error computation."; | |
| 262 | 261 | throw GException::invalid_value(G_RUN, msg); |
| 263 | 262 | } |
| 264 | 263 | |
| 265 | - int npars = m_skymodel->spectral()->size(); | |
| 264 | + // Get number of spectral parameters | |
| 265 | + int npars = model->spectral()->size(); | |
| 266 | 266 | |
| 267 | + // Loop over spectral parameters of sky model | |
| 268 | + for (int i = 0; i < npars; ++i) { | |
| 267 | 269 | |
| 268 | - for (int i = 0; i < npars; ++i) | |
| 269 | - { | |
| 270 | + // Skip parameter if it is fixed | |
| 271 | + if (model->spectral()->at(i).is_fixed()) { | |
| 272 | + continue; | |
| 273 | + } | |
| 270 | 274 | |
| 271 | - // Re-optimize | |
| 272 | - m_obs.optimize(*m_opt); | |
| 275 | + // Initialise with best fitting models | |
| 276 | + m_obs.models(models_best); | |
| 273 | 277 | |
| 274 | - if(m_skymodel->spectral()->at(i).is_fixed()) | |
| 275 | - continue; | |
| 278 | + // Get pointer on sky model parameter | |
| 279 | + GModels& current_models = const_cast<GModels&>(m_obs.models()); | |
| 280 | + GModelSky* current_model = dynamic_cast<GModelSky*>(current_models[m_srcname]); | |
| 281 | + m_model_par = &(current_model->spectral()->at(i)); | |
| 276 | 282 | |
| 277 | - m_model_par = &(m_skymodel->spectral()->at(i)); | |
| 283 | + // Extract current value | |
| 284 | + m_value = m_model_par->factor_value(); | |
| 278 | 285 | |
| 279 | 286 | // Compute parameter bracketing |
| 280 | - m_value = m_model_par->factor_value(); | |
| 281 | 287 | double parmin = std::max(m_model_par->factor_min(), |
| 282 | 288 | m_value - 10*m_model_par->factor_error()); |
| 283 | 289 | double parmax = std::min(m_model_par->factor_max(), |
| ... | ... | @@ -288,7 +294,7 @@ void cterror::run(void) |
| 288 | 294 | log << std::endl; |
| 289 | 295 | log.header1("Compute parameter error"); |
| 290 | 296 | log << gammalib::parformat("Model name"); |
| 291 | - log << m_skymodel->name() << std::endl; | |
| 297 | + log << model->name() << std::endl; | |
| 292 | 298 | log << gammalib::parformat("Parameter name"); |
| 293 | 299 | log << m_model_par->name() << std::endl; |
| 294 | 300 | log << gammalib::parformat("Confidence level"); |
| ... | ... | @@ -304,43 +310,38 @@ void cterror::run(void) |
| 304 | 310 | } |
| 305 | 311 | |
| 306 | 312 | // Compute lower boundary |
| 307 | - error_bisection(parmin, m_value); | |
| 308 | - double value_lo = m_model_par->factor_value(); | |
| 313 | + double value_lo = error_bisection(parmin, m_value); | |
| 309 | 314 | |
| 315 | + // Dump lower parameter value | |
| 310 | 316 | if (logTerse()) { |
| 311 | 317 | log << gammalib::parformat("Final lower parameter value"); |
| 312 | 318 | log << value_lo << std::endl; |
| 313 | 319 | } |
| 314 | 320 | |
| 315 | - // Reset parameter to original value | |
| 316 | - m_model_par->factor_value(m_value); | |
| 317 | - | |
| 318 | 321 | // Compute upper boundary |
| 319 | - error_bisection(m_value, parmax); | |
| 320 | - double value_hi = m_model_par->factor_value(); | |
| 322 | + double value_hi = error_bisection(m_value, parmax); | |
| 321 | 323 | |
| 324 | + // Dump upper parameter value | |
| 322 | 325 | if (logTerse()) { |
| 323 | 326 | log << gammalib::parformat("Final upper parameter"); |
| 324 | 327 | log << value_hi << std::endl; |
| 325 | 328 | } |
| 326 | 329 | |
| 327 | - // Reset parameter to original value | |
| 328 | - m_model_par->factor_value(m_value); | |
| 329 | - | |
| 330 | 330 | // Compute errors |
| 331 | - m_error = 0.5*(value_hi - value_lo); | |
| 332 | - double error_max = std::max(value_hi-m_value, m_value-value_lo); | |
| 333 | - double error_min = std::min(value_hi-m_value, m_value-value_lo); | |
| 334 | - | |
| 331 | + m_error = 0.5 * (value_hi - value_lo); | |
| 332 | + double error_neg = value_lo - m_value; | |
| 333 | + double error_pos = value_hi - m_value; | |
| 334 | + //double error_max = std::max(value_hi-m_value, m_value-value_lo); | |
| 335 | + //double error_min = std::min(value_hi-m_value, m_value-value_lo); | |
| 335 | 336 | |
| 336 | 337 | // Write results to logfile |
| 337 | 338 | if (logTerse()) { |
| 338 | 339 | log << std::endl; |
| 339 | 340 | log.header1("Error results"); |
| 340 | - log << gammalib::parformat("Max error"); | |
| 341 | - log << error_max; | |
| 342 | - log << gammalib::parformat("Min error"); | |
| 343 | - log << error_min; | |
| 341 | + log << gammalib::parformat("Negative error"); | |
| 342 | + log << error_neg; | |
| 343 | + log << gammalib::parformat("Positive error"); | |
| 344 | + log << error_pos; | |
| 344 | 345 | log << gammalib::parformat("Mean parameter error"); |
| 345 | 346 | log << m_error; |
| 346 | 347 | log << " " << m_model_par->unit(); |
| ... | ... | @@ -351,24 +352,24 @@ void cterror::run(void) |
| 351 | 352 | } |
| 352 | 353 | |
| 353 | 354 | // Save error result |
| 354 | - m_model_par->factor_error(m_error); | |
| 355 | + model->spectral()->at(i).factor_error(m_error); | |
| 355 | 356 | |
| 356 | - } | |
| 357 | + } // endfor: looped over spectral parameters | |
| 357 | 358 | |
| 358 | - } | |
| 359 | + // Restore best fitting models (now with new errors computed) | |
| 360 | + m_obs.models(models_best); | |
| 359 | 361 | |
| 362 | + } // endif: source model exists | |
| 360 | 363 | |
| 361 | 364 | // Return |
| 362 | 365 | return; |
| 363 | 366 | } |
| 364 | 367 | |
| 368 | + | |
| 365 | 369 | /***********************************************************************//** |
| 366 | - * @brief Save maps | |
| 370 | + * @brief Save source model | |
| 367 | 371 | * |
| 368 | - * This method saves the error to an ascii file. | |
| 369 | - * | |
| 370 | - * @todo No yet implemented as we have no clear use case yet for saving | |
| 371 | - * the result. | |
| 372 | + * This method saves the soutce model to an XML file. | |
| 372 | 373 | ***************************************************************************/ |
| 373 | 374 | void cterror::save(void) |
| 374 | 375 | { |
| ... | ... | @@ -378,9 +379,8 @@ void cterror::save(void) |
| 378 | 379 | log.header1("Save results"); |
| 379 | 380 | } |
| 380 | 381 | |
| 381 | - | |
| 382 | 382 | // Get output filename |
| 383 | - std::string outmodel = (*this)["outfile"].filename(); | |
| 383 | + std::string outmodel = (*this)["outmodel"].filename(); | |
| 384 | 384 | |
| 385 | 385 | // Write results out as XML model |
| 386 | 386 | if (gammalib::toupper(outmodel) != "NONE") { |
| ... | ... | @@ -391,6 +391,7 @@ void cterror::save(void) |
| 391 | 391 | return; |
| 392 | 392 | } |
| 393 | 393 | |
| 394 | + | |
| 394 | 395 | /*========================================================================== |
| 395 | 396 | = = |
| 396 | 397 | = Private methods = |
| ... | ... | @@ -405,8 +406,6 @@ void cterror::init_members(void) |
| 405 | 406 | // Initialise user parameters |
| 406 | 407 | m_srcname.clear(); |
| 407 | 408 | m_confidence = 0.68; |
| 408 | - m_sigma_min = 0.0; | |
| 409 | - m_sigma_max = 0.0; | |
| 410 | 409 | m_tol = 1.0e-3; |
| 411 | 410 | m_max_iter = 50; |
| 412 | 411 | m_value = 0.0; |
| ... | ... | @@ -415,9 +414,8 @@ void cterror::init_members(void) |
| 415 | 414 | // Initialise protected members |
| 416 | 415 | m_obs.clear(); |
| 417 | 416 | m_dlogL = 0.0; |
| 418 | - m_skymodel = NULL; | |
| 419 | - m_model_par = NULL; | |
| 420 | 417 | m_best_logL = 0.0; |
| 418 | + m_model_par = NULL; | |
| 421 | 419 | |
| 422 | 420 | // Return |
| 423 | 421 | return; |
| ... | ... | @@ -434,15 +432,14 @@ void cterror::copy_members(const cterror& app) |
| 434 | 432 | // Copy user parameters |
| 435 | 433 | m_srcname = app.m_srcname; |
| 436 | 434 | m_confidence = app.m_confidence; |
| 437 | - m_sigma_min = app.m_sigma_min; | |
| 438 | - m_sigma_max = app.m_sigma_max; | |
| 439 | 435 | m_tol = app.m_tol; |
| 440 | 436 | m_max_iter = app.m_max_iter; |
| 441 | 437 | |
| 442 | 438 | // Copy protected members |
| 443 | - m_obs = app.m_obs; | |
| 444 | - m_dlogL = app.m_dlogL; | |
| 445 | - m_best_logL = app.m_best_logL; | |
| 439 | + m_obs = app.m_obs; | |
| 440 | + m_dlogL = app.m_dlogL; | |
| 441 | + m_best_logL = app.m_best_logL; | |
| 442 | + m_model_par = NULL; | |
| 446 | 443 | |
| 447 | 444 | // Return |
| 448 | 445 | return; |
| ... | ... | @@ -506,19 +503,13 @@ void cterror::get_parameters(void) |
| 506 | 503 | |
| 507 | 504 | // Get confidence level and transform into log-likelihood difference |
| 508 | 505 | m_confidence = (*this)["confidence"].real(); |
| 509 | - | |
| 510 | 506 | double sigma = gammalib::erfinv(m_confidence) * gammalib::sqrt_two; |
| 511 | 507 | m_dlogL = (sigma*sigma) / 2.0; |
| 512 | 508 | |
| 513 | - // Read starting boundaries for bisection | |
| 514 | - m_sigma_min = (*this)["sigma_min"].real(); | |
| 515 | - m_sigma_max = (*this)["sigma_max"].real(); | |
| 516 | - | |
| 517 | - // Read precision | |
| 509 | + // Read other parameters | |
| 518 | 510 | m_tol = (*this)["tol"].real(); |
| 519 | 511 | m_max_iter = (*this)["max_iter"].integer(); |
| 520 | 512 | |
| 521 | - | |
| 522 | 513 | // Return |
| 523 | 514 | return; |
| 524 | 515 | } |
| ... | ... | @@ -532,7 +523,7 @@ void cterror::get_parameters(void) |
| 532 | 523 | * |
| 533 | 524 | * This method calculates the error using a bisection method. |
| 534 | 525 | ***************************************************************************/ |
| 535 | -void cterror::error_bisection(const double& min, const double& max) | |
| 526 | +double cterror::error_bisection(const double& min, const double& max) | |
| 536 | 527 | { |
| 537 | 528 | // Copy values to working values |
| 538 | 529 | double wrk_min = min; |
| ... | ... | @@ -541,6 +532,9 @@ void cterror::error_bisection(const double& min, const double& max) |
| 541 | 532 | // Initialise iteration counter |
| 542 | 533 | int iter = 0; |
| 543 | 534 | |
| 535 | + // Initialize mid value | |
| 536 | + double mid = (wrk_min + wrk_max) / 2.0; | |
| 537 | + | |
| 544 | 538 | // Loop until breaking condition is reached |
| 545 | 539 | while (true) { |
| 546 | 540 | |
| ... | ... | @@ -556,8 +550,7 @@ void cterror::error_bisection(const double& min, const double& max) |
| 556 | 550 | |
| 557 | 551 | // Throw exception if maximum iterations are reached |
| 558 | 552 | if (iter > m_max_iter) { |
| 559 | - if(wrk_min - m_model_par->factor_min() < m_tol) | |
| 560 | - { | |
| 553 | + if(wrk_min - m_model_par->factor_min() < m_tol) { | |
| 561 | 554 | std::string msg = "The "+m_model_par->name()+" parameter minimum has been" |
| 562 | 555 | " reached during error calculation. To obtain accurate " |
| 563 | 556 | " errors, consider to set the minimum to a lower value," |
| ... | ... | @@ -567,8 +560,7 @@ void cterror::error_bisection(const double& min, const double& max) |
| 567 | 560 | } |
| 568 | 561 | break; |
| 569 | 562 | } |
| 570 | - else if(m_model_par->factor_max() - wrk_max < m_tol) | |
| 571 | - { | |
| 563 | + else if(m_model_par->factor_max() - wrk_max < m_tol) { | |
| 572 | 564 | std::string msg = "The "+m_model_par->name()+" parameter maximum has been" |
| 573 | 565 | " reached during error calculation. To obtain accurate " |
| 574 | 566 | " errors, consider to set the maxmimum to a higher value" |
| ... | ... | @@ -578,8 +570,7 @@ void cterror::error_bisection(const double& min, const double& max) |
| 578 | 570 | } |
| 579 | 571 | break; |
| 580 | 572 | } |
| 581 | - else | |
| 582 | - { | |
| 573 | + else { | |
| 583 | 574 | std::string msg = "The maximum number of "+gammalib::str(m_max_iter)+ |
| 584 | 575 | " has been reached. You may consider to increase" |
| 585 | 576 | " the \"max_iter\" parameter, and re-run cterror."; |
| ... | ... | @@ -588,7 +579,7 @@ void cterror::error_bisection(const double& min, const double& max) |
| 588 | 579 | } |
| 589 | 580 | |
| 590 | 581 | // Compute center of boundary |
| 591 | - double mid = (wrk_min + wrk_max) / 2.0; | |
| 582 | + mid = (wrk_min + wrk_max) / 2.0; | |
| 592 | 583 | |
| 593 | 584 | // Calculate function value |
| 594 | 585 | double eval_mid = evaluate(mid); |
| ... | ... | @@ -598,9 +589,9 @@ void cterror::error_bisection(const double& min, const double& max) |
| 598 | 589 | break; |
| 599 | 590 | } |
| 600 | 591 | |
| 601 | - // If we are on the crescent side of the parabola | |
| 602 | - if (mid > m_value) | |
| 603 | - { | |
| 592 | + // If we are on the crescent side of the parabola ... | |
| 593 | + if (mid > m_value) { | |
| 594 | + | |
| 604 | 595 | // Change boundaries for further iteration |
| 605 | 596 | if (eval_mid > 0.0) { |
| 606 | 597 | wrk_max = mid; |
| ... | ... | @@ -609,9 +600,10 @@ void cterror::error_bisection(const double& min, const double& max) |
| 609 | 600 | wrk_min = mid; |
| 610 | 601 | } |
| 611 | 602 | } |
| 612 | - // If we are on the decrescent side of the parabola | |
| 613 | - else | |
| 614 | - { | |
| 603 | + | |
| 604 | + // ... otherwise we are on the decrescent side of the parabola | |
| 605 | + else { | |
| 606 | + | |
| 615 | 607 | // Change boundaries for further iteration |
| 616 | 608 | if (eval_mid > 0.0) { |
| 617 | 609 | wrk_min = mid; |
| ... | ... | @@ -626,8 +618,8 @@ void cterror::error_bisection(const double& min, const double& max) |
| 626 | 618 | |
| 627 | 619 | } // endwhile |
| 628 | 620 | |
| 629 | - // Return | |
| 630 | - return; | |
| 621 | + // Return mid value | |
| 622 | + return mid; | |
| 631 | 623 | |
| 632 | 624 | } |
| 633 | 625 | |
| ... | ... | @@ -652,15 +644,12 @@ double cterror::evaluate(const double& value) |
| 652 | 644 | // Change parameter factor |
| 653 | 645 | m_model_par->factor_value(value); |
| 654 | 646 | |
| 655 | - //Fix parameter | |
| 647 | + // Fix parameter | |
| 656 | 648 | m_model_par->fix(); |
| 657 | 649 | |
| 658 | 650 | // Re-optimize |
| 659 | 651 | m_obs.optimize(*m_opt); |
| 660 | 652 | |
| 661 | - // Evaluate likelihood for new model container | |
| 662 | - m_obs.eval(); | |
| 663 | - | |
| 664 | 653 | // Retrieve likelihood |
| 665 | 654 | logL = m_obs.logL(); |
| 666 | 655 | |
| ... | ... | @@ -673,7 +662,7 @@ double cterror::evaluate(const double& value) |
| 673 | 662 | else { |
| 674 | 663 | std::string msg = "Value of parameter \""+m_model_par->name()+"\" " |
| 675 | 664 | "outside of validity range requested. To omit " |
| 676 | - "this error please enlarge the parameter range " | |
| 665 | + "this error please enlarge the parameter range " | |
| 677 | 666 | "in the model XML file."; |
| 678 | 667 | throw GException::invalid_value(G_EVALUATE, msg); |
| 679 | 668 | } |
| ... | ... | @@ -681,7 +670,6 @@ double cterror::evaluate(const double& value) |
| 681 | 670 | // Compute function value |
| 682 | 671 | double logL_difference = logL - m_best_logL - m_dlogL; |
| 683 | 672 | |
| 684 | - | |
| 685 | 673 | // Return |
| 686 | 674 | return logL_difference; |
| 687 | 675 | } | ... | ... |
src/cterror/cterror.hpp
| 1 | 1 | /*************************************************************************** |
| 2 | - * cterror - Parameter error calculation tool * | |
| 2 | + * cterror - Parameter error calculation tool * | |
| 3 | 3 | * ----------------------------------------------------------------------- * |
| 4 | 4 | * copyright (C) 2015 by Florent Forest * |
| 5 | 5 | * ----------------------------------------------------------------------- * |
| ... | ... | @@ -70,13 +70,11 @@ protected: |
| 70 | 70 | void free_members(void); |
| 71 | 71 | void get_parameters(void); |
| 72 | 72 | double evaluate(const double& value); |
| 73 | - void error_bisection(const double& min, const double& max); | |
| 73 | + double error_bisection(const double& min, const double& max); | |
| 74 | 74 | |
| 75 | 75 | // User parameters |
| 76 | 76 | std::string m_srcname; //!< Name of source which is moved around |
| 77 | 77 | double m_confidence; //!< Confidence level |
| 78 | - double m_sigma_min; //!< Starting value minimum (multiple fit errors above fit values) | |
| 79 | - double m_sigma_max; //!< Starting value maximum (multiple fit errors above fit values) | |
| 80 | 78 | double m_tol; //!< Tolerance for limit determination |
| 81 | 79 | int m_max_iter; //!< Maximum number of iterations |
| 82 | 80 | double m_value; //!< Parameter value |
| ... | ... | @@ -85,7 +83,6 @@ protected: |
| 85 | 83 | // Protected members |
| 86 | 84 | GObservations m_obs; //!< Observation container |
| 87 | 85 | double m_dlogL; //!< Likelihood difference for upper limit computation |
| 88 | - GModelSky* m_skymodel; //!< Pointer to sky model | |
| 89 | 86 | GModelPar* m_model_par; //!< Pointer to model parameter |
| 90 | 87 | double m_best_logL; //!< Best fit log likelihood of given model |
| 91 | 88 | GOptimizerLM* m_opt; | ... | ... |
src/cterror/cterror.par
| ... | ... | @@ -19,20 +19,18 @@ |
| 19 | 19 | #=================== |
| 20 | 20 | inobs, f, a, "events.fits",,, "Input event list, counts cube or observation definition file" |
| 21 | 21 | inmodel, f, a, "$CTOOLS/share/models/crab.xml",,, "Source model" |
| 22 | -outfile, f, a, "results.xml",,, "Save results" | |
| 22 | +outmodel, f, a, "results.xml",,, "Source model output file" | |
| 23 | 23 | srcname, s, a, "Crab",,, "Source of interest" |
| 24 | 24 | expcube, f, a, "NONE",,, "Exposure cube file (only needed for stacked analysis)" |
| 25 | 25 | psfcube, f, a, "NONE",,, "PSF cube file (only needed for stacked analysis)" |
| 26 | 26 | bkgcube, f, a, "NONE",,, "Background cube file (only needed for stacked analysis)" |
| 27 | -caldb, s, a, "dummy",,, "Calibration database" | |
| 28 | -irf, s, a, "cta_dummy_irf",,, "Instrument response function" | |
| 27 | +caldb, s, a, "prod2",,, "Calibration database" | |
| 28 | +irf, s, a, "South_50h",,, "Instrument response function" | |
| 29 | 29 | |
| 30 | 30 | # |
| 31 | 31 | # Upper limit calculation parameters |
| 32 | 32 | #=================================== |
| 33 | 33 | confidence, r, h, 0.68,0.0,1.0, "Confidence level" |
| 34 | -sigma_min, r, h, 0.0,,, "Start minimum value (multiple sigmas above fit value)" | |
| 35 | -sigma_max, r, h, 10.0,,, "Start maximum value (multiple sigmas above fit value)" | |
| 36 | 34 | tol, r, h, 1e-3,,, "Computation tolerance" |
| 37 | 35 | max_iter, i, h, 50,,, "Maximum number of iterations" |
| 38 | 36 | ... | ... |
test/test_cscripts.sh.in
| ... | ... | @@ -42,6 +42,7 @@ rm -rf *.fits *.xml *.dat |
| 42 | 42 | # |
| 43 | 43 | # Creates pfiles directory |
| 44 | 44 | # ======================== |
| 45 | +rm -rf pfiles | |
| 45 | 46 | mkdir -p pfiles |
| 46 | 47 | export PFILES=pfiles |
| 47 | 48 | |
| ... | ... | @@ -316,3 +317,5 @@ else |
| 316 | 317 | exit 1 |
| 317 | 318 | fi |
| 318 | 319 | $ECHO " ok" |
| 320 | + | |
| 321 | +ls -la pfiles | ... | ... |
test/test_ctools.sh
| ... | ... | @@ -47,6 +47,7 @@ ctbkgcube=../src/ctbkgcube/ctbkgcube |
| 47 | 47 | ctcubemask=../src/ctcubemask/ctcubemask |
| 48 | 48 | ctbutterfly=../src/ctbutterfly/ctbutterfly |
| 49 | 49 | ctulimit=../src/ctulimit/ctulimit |
| 50 | +cterror=../src/cterror/cterror | |
| 50 | 51 | |
| 51 | 52 | |
| 52 | 53 | # |
| ... | ... | @@ -61,7 +62,7 @@ rm -rf ulimit.dat butterfly.txt |
| 61 | 62 | # ======================== |
| 62 | 63 | rm -rf pfiles |
| 63 | 64 | mkdir -p pfiles |
| 64 | -cp -r ../src/*/*.par pfiles/ | |
| 65 | +#cp -r ../src/*/*.par pfiles/ | |
| 65 | 66 | export PFILES=pfiles |
| 66 | 67 | |
| 67 | 68 | |
| ... | ... | @@ -595,3 +596,24 @@ else |
| 595 | 596 | fi |
| 596 | 597 | $ECHO " ok" |
| 597 | 598 | |
| 599 | + | |
| 600 | +# | |
| 601 | +# Test cterror | |
| 602 | +# ============ | |
| 603 | +$ECHO -n "Test cterror: " | |
| 604 | +$cterror inobs="data/crab_events.fits.gz" \ | |
| 605 | + inmodel="data/crab.xml" \ | |
| 606 | + outmodel="results_error.xml" \ | |
| 607 | + srcname="Crab" \ | |
| 608 | + caldb="irf" \ | |
| 609 | + irf="cta_dummy_irf" | |
| 610 | +$ECHO -n "." | |
| 611 | +if [ -s "cterror.log" ] | |
| 612 | +then | |
| 613 | + $ECHO -n "." | |
| 614 | +else | |
| 615 | + $ECHO " cterror.log file is not found" | |
| 616 | + exit 1 | |
| 617 | +fi | |
| 618 | +$ECHO " ok" | |
| 619 | + | ... | ... |
test/test_python.py
| ... | ... | @@ -2,7 +2,7 @@ |
| 2 | 2 | # ========================================================================== |
| 3 | 3 | # This scripts performs unit tests for the ctools package. |
| 4 | 4 | # |
| 5 | -# Copyright (C) 2012-2014 Juergen Knoedlseder | |
| 5 | +# Copyright (C) 2012-2015 Juergen Knoedlseder | |
| 6 | 6 | # |
| 7 | 7 | # This program is free software: you can redistribute it and/or modify |
| 8 | 8 | # it under the terms of the GNU General Public License as published by |
| ... | ... | @@ -35,6 +35,7 @@ import test_ctbkgcube |
| 35 | 35 | import test_ctcubemask |
| 36 | 36 | import test_ctbutterfly |
| 37 | 37 | import test_ctulimit |
| 38 | +import test_cterror | |
| 38 | 39 | import test_pipelines |
| 39 | 40 | |
| 40 | 41 | |
| ... | ... | @@ -62,6 +63,7 @@ if __name__ == '__main__': |
| 62 | 63 | suite_ctcubemask = test_ctcubemask.Test() |
| 63 | 64 | suite_ctbutterfly = test_ctbutterfly.Test() |
| 64 | 65 | suite_ctulimit = test_ctulimit.Test() |
| 66 | + suite_cterror = test_cterror.Test() | |
| 65 | 67 | suite_pipelines = test_pipelines.Test() |
| 66 | 68 | |
| 67 | 69 | # Setup unit tests |
| ... | ... | @@ -78,6 +80,7 @@ if __name__ == '__main__': |
| 78 | 80 | suite_ctcubemask.set() |
| 79 | 81 | suite_ctbutterfly.set() |
| 80 | 82 | suite_ctulimit.set() |
| 83 | + suite_cterror.set() | |
| 81 | 84 | suite_pipelines.set() |
| 82 | 85 | |
| 83 | 86 | # Append tests to container |
| ... | ... | @@ -94,6 +97,7 @@ if __name__ == '__main__': |
| 94 | 97 | suites.append(suite_ctcubemask) |
| 95 | 98 | suites.append(suite_ctbutterfly) |
| 96 | 99 | suites.append(suite_ctulimit) |
| 100 | + suites.append(suite_cterror) | |
| 97 | 101 | suites.append(suite_pipelines) |
| 98 | 102 | |
| 99 | 103 | # Set PFILES environment variable | ... | ... |