Michael Mayer / ctools

Commit 79c614e4061244739267731d313d6208a3bc019b

Authored by Jürgen Knödlseder
2 parents 06535a17 edef6971

Merge branch '1832-ulimit-new-plaw-parameters' into integration

Inline Side-by-side
Showing 2 changed files with 113 additions and 152 deletions
src/ctulimit/ctulimit.cpp
  View file @79c614e
... ... @@ -241,17 +241,8 @@ void ctulimit::run(void)
241 241 }
242 242 }
243 243  
244   - // Write observation(s) into logger
245   - if (logTerse()) {
246   - log << std::endl;
247   - if (m_obs.size() > 1) {
248   - log.header1("Observations");
249   - }
250   - else {
251   - log.header1("Observation");
252   - }
253   - log << m_obs << std::endl;
254   - }
  244 + // Write input observation container into logger
  245 + log_observations(NORMAL, m_obs, "Input observation");
255 246  
256 247 // Save original models
257 248 GModels models_orig = m_obs.models();
... ... @@ -264,27 +255,21 @@ void ctulimit::run(void)
264 255 m_best_logL = m_obs.logL();
265 256 if (m_best_logL == 0.0) {
266 257  
267   - // Write header
268   - if (logTerse()) {
269   - log << std::endl;
270   - log.header1("Compute best-fit likelihood");
271   - }
  258 + // Write header into logger
  259 + log_header1(TERSE, "Compute best-fit likelihood");
272 260  
273   - // Optimize and save best log-likelihood
274   - m_obs.optimize(m_opt);
275   - m_obs.errors(m_opt);
276   - m_best_logL = m_obs.logL();
  261 + // Optimize and save best log-likelihood
  262 + m_obs.optimize(m_opt);
  263 + m_obs.errors(m_opt);
  264 + m_best_logL = m_obs.logL();
277 265  
278   - // Store optimizer for later recovery
279   - best_opt = m_opt;
  266 + // Store optimizer for later recovery
  267 + best_opt = m_opt;
280 268  
281   - // Write optimised model into logger
282   - if (logTerse()) {
283   - log << m_opt << std::endl;
284   - log << gammalib::parformat("Maximum log likelihood");
285   - log << gammalib::str(m_best_logL,3) << std::endl;
286   - log << m_obs.models() << std::endl;
287   - }
  269 + // Write optimised model into logger
  270 + log_string(NORMAL, m_opt.print(m_chatter));
  271 + log_value(NORMAL,"Maximum log likelihood",gammalib::str(m_best_logL,3));
  272 + log_string(NORMAL, m_obs.models().print(m_chatter));
288 273  
289 274 } // endif: likelihood was zero
290 275  
... ... @@ -303,34 +288,22 @@ void ctulimit::run(void)
303 288 double parmax = std::min(m_model_par->factor_max(),
304 289 value + m_sigma_max * error);
305 290  
306   - // Write header
307   - if (logTerse()) {
308   - log << std::endl;
309   - log.header1("Compute upper limit");
310   - log << gammalib::parformat("Model name");
311   - log << m_skymodel->name() << std::endl;
312   - log << gammalib::parformat("Parameter name");
313   - log << m_model_par->name() << std::endl;
314   - log << gammalib::parformat("Confidence level");
315   - log << m_confidence*100.0 << "%" << std::endl;
316   - log << gammalib::parformat("Log-likelihood difference");
317   - log << m_dlogL << std::endl;
318   - log << gammalib::parformat("Initial parameter range");
319   - log << "[";
320   - log << parmin;
321   - log << ", ";
322   - log << parmax;
323   - log << "]" << std::endl;
324   - }
  291 + // Write header into logger
  292 + log_header1(TERSE, "Compute upper limit");
  293 +
  294 + // Write some parameters into logger
  295 + log_value(NORMAL, "Model name", m_skymodel->name());
  296 + log_value(NORMAL, "Parameter name", m_model_par->name());
  297 + log_value(NORMAL, "Confidence level", gammalib::str(m_confidence*100.0)+"%");
  298 + log_value(NORMAL, "Log-likelihood difference", m_dlogL);
  299 + log_value(NORMAL, "Initial parameter range",
  300 + "["+gammalib::str(parmin)+", "+gammalib::str(parmax)+"]");
325 301  
326 302 // Compute upper limit
327 303 ulimit_bisection(parmin, parmax);
328 304  
329   - // Write final parameter range
330   - if (logTerse()) {
331   - log << gammalib::parformat("Final parameter");
332   - log << m_model_par->value() << std::endl;
333   - }
  305 + // Write final parameter into logger
  306 + log_value(NORMAL, "Final parameter", m_model_par->value());
334 307  
335 308 // Get reference energy for differential upper limit
336 309 GEnergy eref = GEnergy(m_eref, "TeV");
... ... @@ -344,28 +317,21 @@ void ctulimit::run(void)
344 317 m_flux_ulimit = m_skymodel->spectral()->flux(emin, emax);
345 318 m_eflux_ulimit = m_skymodel->spectral()->eflux(emin, emax);
346 319  
347   - // Write results to logfile
348   - if (logTerse()) {
349   - log << std::endl;
350   - log.header1("Upper limit results");
351   - log << gammalib::parformat("Differential flux limit");
352   - log << m_diff_ulimit;
353   - log << " ph/cm2/s/MeV at ";
354   - log << m_eref << " TeV";
355   - log << std::endl;
356   - log << gammalib::parformat("Integral flux limit");
357   - log << m_flux_ulimit;
358   - log << " ph/cm2/s within [";
359   - log << m_emin << "-" << m_emax;
360   - log << "] TeV";
361   - log << std::endl;
362   - log << gammalib::parformat("Energy flux limit");
363   - log << m_eflux_ulimit;
364   - log << " erg/cm2/s within [";
365   - log << m_emin << "-" << m_emax;
366   - log << "] TeV";
367   - log << std::endl;
368   - }
  320 + // Write header into logger
  321 + log_header1(TERSE, "Upper limit results");
  322 +
  323 + // Write result into logger
  324 + log_value(TERSE, "Differential flux limit",
  325 + gammalib::str(m_diff_ulimit)+" ph/cm2/s/MeV at "+
  326 + gammalib::str(m_eref)+" TeV");
  327 + log_value(TERSE, "Integral flux limit",
  328 + gammalib::str(m_flux_ulimit)+" ph/cm2/s within ["+
  329 + gammalib::str(m_emin)+"-"+
  330 + gammalib::str(m_emax)+"] TeV");
  331 + log_value(TERSE, "Energy flux limit",
  332 + gammalib::str(m_eflux_ulimit)+" erg/cm2/s within ["+
  333 + gammalib::str(m_emin)+"-"+
  334 + gammalib::str(m_emax)+"] TeV");
369 335  
370 336 // Recover original models
371 337 m_obs.models(models_orig);
... ... @@ -396,27 +362,6 @@ void ctulimit::run(void)
396 362 ***************************************************************************/
397 363 void ctulimit::save(void)
398 364 {
399   - /*
400   - // Write header
401   - if (logTerse()) {
402   - log << std::endl;
403   - log.header1("Save upper limit");
404   - }
405   -
406   - // Get output filename
407   - m_outfile = (*this)["outfile"].filename();
408   -
409   - // Create CSV table with 3 columns
410   - GCsv table(1, 3);
411   -
412   - table.real(0, 0, m_diff_ulimit);
413   - table.real(0, 1, m_flux_ulimit);
414   - table.real(0, 2, m_eflux_ulimit);
415   -
416   - // Save CSV table
417   - table.save(m_outfile, " ", clobber());
418   - */
419   -
420 365 // Return
421 366 return;
422 367 }
... ... @@ -444,6 +389,7 @@ void ctulimit::init_members(void)
444 389 m_tol = 1.0e-6;
445 390 m_max_iter = 50;
446 391 m_apply_edisp = false;
  392 + m_chatter = static_cast<GChatter>(2);
447 393  
448 394 // Initialise protected members
449 395 m_obs.clear();
... ... @@ -483,6 +429,7 @@ void ctulimit::copy_members(const ctulimit&amp; app)
483 429 m_tol = app.m_tol;
484 430 m_max_iter = app.m_max_iter;
485 431 m_apply_edisp = app.m_apply_edisp;
  432 + m_chatter = app.m_chatter;
486 433  
487 434 // Copy protected members
488 435 m_obs = app.m_obs;
... ... @@ -547,14 +494,9 @@ void ctulimit::get_parameters(void)
547 494  
548 495 } // endif: no models were associated with observations
549 496  
550   - // Get name of test source and check container for this name
  497 + // Get name of test source. The get_model_parameter() method will check
  498 + // where a source with the given name exists in the model container.
551 499 m_srcname = (*this)["srcname"].string();
552   - if (!m_obs.models().contains(m_srcname)) {
553   - std::string msg = "Source \""+m_srcname+"\" not found in model "
554   - "container. Please add a source with that name "
555   - "or check for possible typos.";
556   - throw GException::invalid_value(G_GET_PARAMETERS, msg);
557   - }
558 500  
559 501 // Get relevant model and parameter for upper limit computation
560 502 get_model_parameter();
... ... @@ -580,6 +522,9 @@ void ctulimit::get_parameters(void)
580 522 m_tol = (*this)["tol"].real();
581 523 m_max_iter = (*this)["max_iter"].integer();
582 524  
  525 + // Get remaining parameters
  526 + m_chatter = static_cast<GChatter>((*this)["chatter"].integer());
  527 +
583 528 // Write parameters into logger
584 529 log_parameters(TERSE);
585 530  
... ... @@ -592,52 +537,73 @@ void ctulimit::get_parameters(void)
592 537 * @brief Get application parameters
593 538 *
594 539 * @exception GException::invalid_value
595   - * Did not find a valid model parameter
  540 + * Did not find a valid model or model parameter
596 541 *
597 542 * Extracts a pointer to the sky model (m_skymodel) and a pointer to the
598   - * relevant model parameter (m_model_par) from the model container.
  543 + * relevant model parameter (m_model_par) from the model container. If no
599 544 ***************************************************************************/
600 545 void ctulimit::get_model_parameter(void)
601 546 {
602   - // Continue only if source model exists
603   - if (m_obs.models().contains(m_srcname)) {
604   -
605   - // Get relevant model and parameter for upper limit computation.
606   - GModels& models = const_cast<GModels&>(m_obs.models());
607   - m_skymodel = dynamic_cast<GModelSky*>(models[m_srcname]);
608   - if (m_skymodel == NULL) {
609   - std::string msg = "Source \""+m_srcname+"\" is not a sky model. "
610   - "Please specify the name of a sky model for "
611   - "upper limit computation.";
612   - throw GException::invalid_value(G_GET_MODEL_PARAMETER, msg);
613   - }
614   - if (m_skymodel->spectral()->type() == "NodeFunction") {
615   - std::string msg = "\"NodeFunction\" cannot be used as spectral "
616   - "model for an upper limit computation. "
617   - "Please specify another spectral model.";
618   - throw GException::invalid_value(G_GET_MODEL_PARAMETER, msg);
619   - }
620   - if (m_skymodel->spectral()->has_par("Normalization")) {
621   - m_model_par = &(m_skymodel->spectral()->operator[]("Normalization"));
622   - }
623   - else if (m_skymodel->spectral()->has_par("Prefactor")) {
624   - m_model_par = &(m_skymodel->spectral()->operator[]("Prefactor"));
625   - }
626   - else if (m_skymodel->spectral()->has_par("Integral")) {
627   - m_model_par = &(m_skymodel->spectral()->operator[]("Integral"));
628   - }
629   - else if (m_skymodel->spectral()->has_par("Value")) {
630   - m_model_par = &(m_skymodel->spectral()->operator[]("Value"));
631   - }
632   - else {
633   - std::string msg = "Require spectral parameter \"Normalization\", "
634   - "\"Prefactor\", \"Integral\" or \"Value\" for "
635   - "upper limit computation. The specified source "
636   - "\""+m_srcname+"\" does not have such a parameter.";
637   - throw GException::invalid_value(G_GET_MODEL_PARAMETER, msg);
  547 + // Define list of valid spectral model parameters, terminated with a
  548 + // NULL pointer to signal the end of the list
  549 + const char* pars[] = {"Normalization", "Value", "Prefactor",
  550 + "Integral", "PhotonFlux", "EnergyFlux",
  551 + NULL};
  552 +
  553 + // Initialises sky model and model parameters pointers
  554 + m_skymodel = NULL;
  555 + m_model_par = NULL;
  556 +
  557 + // If the model container does not container the specified source then
  558 + // throw an exception
  559 + if (!m_obs.models().contains(m_srcname)) {
  560 + std::string msg = "Source \""+m_srcname+"\" not found in model "
  561 + "container. Please add a source with that name "
  562 + "or check for possible typos.";
  563 + throw GException::invalid_value(G_GET_MODEL_PARAMETER, msg);
  564 + }
  565 +
  566 + // Get relevant sky model for upper limit computation. If the model is not
  567 + // a sky model or if the model has a "NodeFunction" as spectral component
  568 + // then throw an exception.
  569 + GModels& models = const_cast<GModels&>(m_obs.models());
  570 + m_skymodel = dynamic_cast<GModelSky*>(models[m_srcname]);
  571 + if (m_skymodel == NULL) {
  572 + std::string msg = "Source \""+m_srcname+"\" is not a sky model. Please "
  573 + "specify the name of a sky model for upper limit "
  574 + "computation.";
  575 + throw GException::invalid_value(G_GET_MODEL_PARAMETER, msg);
  576 + }
  577 + if (m_skymodel->spectral()->type() == "NodeFunction") {
  578 + std::string msg = "\"NodeFunction\" cannot be used as spectral model "
  579 + "for an upper limit computation. Please specify "
  580 + "another spectral model.";
  581 + throw GException::invalid_value(G_GET_MODEL_PARAMETER, msg);
  582 + }
  583 +
  584 + // Find appropriate model parameter and store its pointer in the
  585 + // m_model_par member
  586 + for (const char** par = pars; *par != NULL; ++par) {
  587 + if (m_skymodel->spectral()->has_par(*par)) {
  588 + m_model_par = &(m_skymodel->spectral()->operator[](*par));
  589 + break;
638 590 }
  591 + }
639 592  
640   - } // endif: source model existed
  593 + // If no model parameter was found then throw an exception
  594 + if (m_model_par == NULL) {
  595 + std::string msg = "Require one of the following spectral "
  596 + "parameters for upper limit computation: ";
  597 + for (const char** par = pars; *par != NULL; ++par) {
  598 + if (par != pars) {
  599 + msg.append(", ");
  600 + }
  601 + msg.append("\""+std::string(*par)+"\"");
  602 + }
  603 + msg.append(". The specified source \""+m_srcname+"\" does not "
  604 + "have such a parameter.");
  605 + throw GException::invalid_value(G_GET_MODEL_PARAMETER, msg);
  606 + }
641 607  
642 608 // Return
643 609 return;
... ... @@ -665,14 +631,8 @@ void ctulimit::ulimit_bisection(const double&amp; min, const double&amp; max)
665 631 while (true) {
666 632  
667 633 // Log information
668   - if (logExplicit()) {
669   - log << gammalib::parformat("Iteration "+gammalib::str(iter));
670   - log << "[";
671   - log << wrk_min;
672   - log << ", ";
673   - log << wrk_max;
674   - log << "]" << std::endl;
675   - }
  634 + log_value(EXPLICIT, "Iteration "+gammalib::str(iter),
  635 + "["+gammalib::str(wrk_min)+", "+gammalib::str(wrk_max)+"]");
676 636  
677 637 // Throw exception if maximum iterations are reached
678 638 if (iter > m_max_iter) {
... ...
src/ctulimit/ctulimit.hpp
  View file @79c614e
1 1 /***************************************************************************
2 2 * ctulimit - Upper limit calculation tool *
3 3 * ----------------------------------------------------------------------- *
4   - * copyright (C) 2015 by Michael Mayer *
  4 + * copyright (C) 2015-2016 by Michael Mayer *
5 5 * ----------------------------------------------------------------------- *
6 6 * *
7 7 * This program is free software: you can redistribute it and/or modify *
... ... @@ -34,7 +34,7 @@
34 34  
35 35 /* __Definitions _________________________________________________________ */
36 36 #define CTULIMIT_NAME "ctulimit"
37   -#define CTULIMIT_VERSION "01-00-00"
  37 +#define CTULIMIT_VERSION "1.2.0"
38 38  
39 39  
40 40 /***********************************************************************//**
... ... @@ -94,6 +94,7 @@ protected:
94 94 double m_tol; //!< Tolerance for limit determination
95 95 int m_max_iter; //!< Maximum number of iterations
96 96 bool m_apply_edisp; //!< Apply energy dispersion?
  97 + GChatter m_chatter; //!< Chattiness
97 98  
98 99 // Protected members
99 100 GObservations m_obs; //!< Observation container
... ...