diff --git a/src/cterror/cterror.cpp b/src/cterror/cterror.cpp index 5063d8a..f477d99 100644 --- a/src/cterror/cterror.cpp +++ b/src/cterror/cterror.cpp @@ -228,35 +228,15 @@ void cterror::run(void) // Get task parameters get_parameters(); - // Set energy dispersion flag for all CTA observations and save old + // Set energy dispersion flags of all CTA observations and save old // values in save_edisp vector - std::vector save_edisp; - save_edisp.assign(m_obs.size(), false); - for (int i = 0; i < m_obs.size(); ++i) { - GCTAObservation* obs = dynamic_cast(m_obs[i]); - if (obs != NULL) { - save_edisp[i] = obs->response()->apply_edisp(); - obs->response()->apply_edisp(m_apply_edisp); - } - } + std::vector save_edisp = set_edisp(m_obs, m_apply_edisp); - // Write observation(s) into logger - if (logTerse()) { - log << std::endl; - if (m_obs.size() > 1) { - log.header1("Observations"); - } - else { - log.header1("Observation"); - } - log << m_obs << std::endl; - } + // Write input observation container into logger + log_observations(NORMAL, m_obs, "Input observation"); - // Write header - if (logTerse()) { - log << std::endl; - log.header1("Compute best-fit likelihood"); - } + // Write header into logger + log_header1(TERSE, "Compute best-fit likelihood"); // Optimize and save best log-likelihood m_obs.optimize(m_opt); @@ -266,13 +246,10 @@ void cterror::run(void) // Store optimizer for later recovery GOptimizerLM best_opt = m_opt; - // Write optimised model into logger - if (logTerse()) { - log << m_opt << std::endl; - log << gammalib::parformat("Maximum log likelihood"); - log << gammalib::str(m_best_logL,3) << std::endl; - log << m_obs.models() << std::endl; - } + // Write optimisation results and models into logger + log_string(NORMAL, m_opt.print(m_chatter)); + log_value(NORMAL, "Maximum log likelihood", gammalib::str(m_best_logL,3)); + log_string(NORMAL, m_obs.models().print(m_chatter)); // Continue only if source model exists if (m_obs.models().contains(m_srcname)) { @@ -310,63 +287,39 @@ void cterror::run(void) double parmax = std::min(m_model_par->factor_max(), m_value + 10.0*m_model_par->factor_error()); - // Write header - if (logTerse()) { - log << std::endl; - log.header1("Compute error for source \""+m_srcname+"\"" - " parameter \""+m_model_par->name()+"\""); - log << gammalib::parformat("Confidence level"); - log << m_confidence*100.0 << "%" << std::endl; - log << gammalib::parformat("Log-likelihood difference"); - log << m_dlogL << std::endl; - log << gammalib::parformat("Initial factor range"); - log << "["; - log << parmin; - log << ", "; - log << parmax; - log << "]" << std::endl; - } + // Write header and initial parameters into logger + log_header1(TERSE, "Compute error for source \""+m_srcname+"\"" + " parameter \""+m_model_par->name()+"\""); + log_value(NORMAL, "Confidence level", + gammalib::str(m_confidence*100.0)+" %"); + log_value(NORMAL, "Log-likelihood difference", m_dlogL); + log_value(NORMAL, "Initial factor range", + "["+gammalib::str(parmin)+", "+gammalib::str(parmax)+"]"); - // Compute lower boundary + // Compute lower and upper boundaries double value_lo = error_bisection(parmin, m_value); - - // Write lower parameter value - if (logTerse()) { - log << gammalib::parformat("Lower parameter factor"); - log << value_lo << std::endl; - } - - // Compute upper boundary double value_hi = error_bisection(m_value, parmax); - // Write upper parameter value - if (logTerse()) { - log << gammalib::parformat("Upper parameter factor"); - log << value_hi << std::endl; - } - // Compute errors - double error = 0.5 * (value_hi - value_lo); - double error_neg = m_value - value_lo; - double error_pos = value_hi - m_value; - //double error_max = std::max(value_hi-m_value, m_value-value_lo); - //double error_min = std::min(value_hi-m_value, m_value-value_lo); - - // Write errors - if (logTerse()) { - log << gammalib::parformat("Error from curvature"); - log << m_model_par->error(); - log << " " << m_model_par->unit() << std::endl; - log << gammalib::parformat("Error from profile"); - log << std::abs(error*m_model_par->scale()); - log << " " << m_model_par->unit() << std::endl; - log << gammalib::parformat("Negative profile error"); - log << std::abs(error_neg*m_model_par->scale()); - log << " " << m_model_par->unit() << std::endl; - log << gammalib::parformat("Positive profile error"); - log << std::abs(error_pos*m_model_par->scale()); - log << " " << m_model_par->unit() << std::endl; - } + double error = 0.5 * (value_hi - value_lo); + double error_neg = m_value - value_lo; + double error_pos = value_hi - m_value; + double error_value = std::abs(error*m_model_par->scale()); + double error_value_neg = std::abs(error_neg*m_model_par->scale()); + double error_value_pos = std::abs(error_pos*m_model_par->scale()); + + // Write results into logger + std::string unit = " " + m_model_par->unit(); + log_value(NORMAL, "Lower parameter factor", value_lo); + log_value(NORMAL, "Upper parameter factor", value_hi); + log_value(NORMAL, "Error from curvature", + gammalib::str(m_model_par->error()) + unit); + log_value(NORMAL, "Error from profile", + gammalib::str(error_value) + unit); + log_value(NORMAL, "Negative profile error", + gammalib::str(error_value_neg) + unit); + log_value(NORMAL, "Positive profile error", + gammalib::str(error_value_pos) + unit); // Save error result model->at(i).factor_error(error); @@ -381,13 +334,8 @@ void cterror::run(void) // Recover optimizer m_opt = best_opt; - // Restore energy dispersion flag for all CTA observations - for (int i = 0; i < m_obs.size(); ++i) { - GCTAObservation* obs = dynamic_cast(m_obs[i]); - if (obs != NULL) { - obs->response()->apply_edisp(save_edisp[i]); - } - } + // Restore energy dispersion flags of all CTA observations + restore_edisp(m_obs, save_edisp); // Return return; @@ -402,27 +350,26 @@ void cterror::run(void) void cterror::save(void) { // Write header - if (logTerse()) { - log << std::endl; - log.header1("Save results"); - } + log_header1(TERSE, "Save results"); // Get output filename m_outmodel = (*this)["outmodel"].filename(); - // Save only if filename is non-empty + // Save only if filename is not empty and not "NONE" if (!m_outmodel.is_empty() && gammalib::toupper(m_outmodel.url()) != "NONE") { // Log filename - if (logTerse()) { - log << "Save errors into file \""+m_outmodel+"\"."; - log << std::endl; - } + log_value(NORMAL, "Model definition file", m_outmodel.url()); // Write results out as XML model - m_obs.models().save(m_outmodel); + m_obs.models().save(m_outmodel.url()); + + } + // ... otherwise signal that file was not save + else { + log_value(NORMAL, "Model definition file", "NONE"); } // Return @@ -448,6 +395,7 @@ void cterror::init_members(void) m_tol = 1.0e-3; m_max_iter = 50; m_value = 0.0; + m_chatter = static_cast(2); // Initialise protected members m_obs.clear(); @@ -480,6 +428,7 @@ void cterror::copy_members(const cterror& app) m_tol = app.m_tol; m_max_iter = app.m_max_iter; m_apply_edisp = app.m_apply_edisp; + m_chatter = app.m_chatter; // Copy protected members m_obs = app.m_obs; @@ -559,6 +508,7 @@ void cterror::get_parameters(void) // Read other parameters m_tol = (*this)["tol"].real(); m_max_iter = (*this)["max_iter"].integer(); + m_chatter = static_cast((*this)["chatter"].integer()); // Read ahead parameters that are only needed when the tool gets // executed @@ -606,9 +556,7 @@ double cterror::error_bisection(const double& min, const double& max) "consider setting the minimum parameter " "value to a lower value, and re-run " "cterror."; - if (logTerse()) { - log << msg; - } + log_string(TERSE, msg); break; } else if (m_model_par->factor_max() - wrk_max < m_tol) { @@ -618,9 +566,7 @@ double cterror::error_bisection(const double& min, const double& max) "consider setting the maximum parameter " "value to a higher value, and re-run " "cterror."; - if (logTerse()) { - log << msg; - } + log_string(TERSE, msg); break; } else { @@ -639,17 +585,9 @@ double cterror::error_bisection(const double& min, const double& max) // Calculate function value double eval_mid = evaluate(mid); - // Log interval - if (logExplicit()) { - log << gammalib::parformat(" Iteration "+gammalib::str(iter)); - log << "["; - log << wrk_min; - log << ", "; - log << wrk_max; - log << "] "; - log << eval_mid; - log << std::endl; - } + // Write interval into logger + log_value(EXPLICIT, " Iteration "+gammalib::str(iter), + "["+gammalib::str(wrk_min)+", "+gammalib::str(wrk_max)+"]"); // Check for convergence inside tolerance if (std::abs(eval_mid) < m_tol) { diff --git a/src/cterror/cterror.hpp b/src/cterror/cterror.hpp index 8db9fa4..249666b 100644 --- a/src/cterror/cterror.hpp +++ b/src/cterror/cterror.hpp @@ -35,7 +35,7 @@ /* __Definitions _________________________________________________________ */ #define CTERROR_NAME "cterror" -#define CTERROR_VERSION "1.1.0" +#define CTERROR_VERSION "1.2.0" /***********************************************************************//** @@ -80,6 +80,7 @@ protected: int m_max_iter; //!< Maximum number of iterations double m_value; //!< Parameter value bool m_apply_edisp; //!< Apply energy dispersion? + GChatter m_chatter; //!< Chattiness // Protected members GObservations m_obs; //!< Observation container @@ -87,7 +88,6 @@ protected: GModelPar* m_model_par; //!< Pointer to model parameter double m_best_logL; //!< Best fit log likelihood of given model GOptimizerLM m_opt; //!< Optimizer - }; diff --git a/src/ctlike/ctlike.cpp b/src/ctlike/ctlike.cpp index a7c3f89..28733ca 100644 --- a/src/ctlike/ctlike.cpp +++ b/src/ctlike/ctlike.cpp @@ -215,24 +215,12 @@ void ctlike::run(void) // Get parameters get_parameters(); - // Set energy dispersion flag for all CTA observations and save old + // Set energy dispersion flags of all CTA observations and save old // values in save_edisp vector - std::vector save_edisp; - save_edisp.assign(m_obs.size(), false); - for (int i = 0; i < m_obs.size(); ++i) { - GCTAObservation* obs = dynamic_cast(m_obs[i]); - if (obs != NULL) { - save_edisp[i] = obs->response()->apply_edisp(); - obs->response()->apply_edisp(m_apply_edisp); - } - } + std::vector save_edisp = set_edisp(m_obs, m_apply_edisp); - // Write observation(s) into logger - if (logTerse()) { - log << std::endl; - log.header1(gammalib::number("Observation",m_obs.size())); - log << m_obs.print(m_chatter) << std::endl; - } + // Write input observation container into logger + log_observations(NORMAL, m_obs, "Input observation"); // Optimize model parameters using LM optimizer optimize_lm(); @@ -303,29 +291,16 @@ void ctlike::run(void) } // Write results into logger - if (logTerse()) { - log << std::endl; - log.header1("Maximum likelihood optimisation results"); - log << m_opt->print(m_chatter) << std::endl; - log << gammalib::parformat("Maximum log likelihood"); - log << gammalib::str(m_logL,3) << std::endl; - log << gammalib::parformat("Observed events (Nobs)"); - log << gammalib::str(num_events, 3) << std::endl; - log << gammalib::parformat("Predicted events (Npred)"); - log << gammalib::str(npred, 3); - log << " (Nobs - Npred = "; - log << gammalib::str(num_events-npred); - log << ")" << std::endl; - log << m_obs.models().print(m_chatter) << std::endl; - } + log_header1(NORMAL, "Maximum likelihood optimisation results"); + log_string(NORMAL, m_opt->print(m_chatter)); + log_value(NORMAL, "Maximum log likelihood", gammalib::str(m_logL,3)); + log_value(NORMAL, "Observed events (Nobs)", gammalib::str(num_events,3)); + log_value(NORMAL, "Predicted events (Npred)", gammalib::str(npred,3)+ + " (Nobs - Npred = "+gammalib::str(num_events-npred)+")"); + log_string(NORMAL, m_obs.models().print(m_chatter)); - // Restore energy dispersion flag for all CTA observations - for (int i = 0; i < m_obs.size(); ++i) { - GCTAObservation* obs = dynamic_cast(m_obs[i]); - if (obs != NULL) { - obs->response()->apply_edisp(save_edisp[i]); - } - } + // Restore energy dispersion flags of all CTA observations + restore_edisp(m_obs, save_edisp); // Return return; @@ -340,10 +315,7 @@ void ctlike::run(void) void ctlike::save(void) { // Write header - if (logTerse()) { - log << std::endl; - log.header1("Save results"); - } + log_header1(TERSE, "Save results"); // Get output filename m_outmodel = (*this)["outmodel"].filename(); @@ -353,16 +325,18 @@ void ctlike::save(void) gammalib::toupper(m_outmodel.url()) != "NONE") { // Log filename - if (logTerse()) { - log << "Save fitted parameters into file \""+m_outmodel+"\"."; - log << std::endl; - } + log_value(NORMAL, "Model definition file", m_outmodel.url()); // Write results out as XML model m_obs.models().save(m_outmodel); } + // ... otherwise signal that file was not save + else { + log_value(NORMAL, "Model definition file", "NONE"); + } + // Return return; } @@ -434,7 +408,7 @@ void ctlike::get_parameters(void) } // Set optimizer logger - if (logTerse()) { + if (logNormal()) { static_cast(m_opt)->logger(&log); } else { @@ -456,12 +430,9 @@ void ctlike::get_parameters(void) ***************************************************************************/ void ctlike::optimize_lm(void) { - // Write Header for optimization and indent for optimizer logging - if (logTerse()) { - log << std::endl; - log.header1("Maximum likelihood optimisation"); - log.indent(1); - } + // Write header + log_header1(TERSE, "Maximum likelihood optimisation"); + log.indent(1); // Compute number of fitted parameters int nfit = 0; @@ -476,14 +447,9 @@ void ctlike::optimize_lm(void) // Notify if all parameters are fixed if (nfit == 0) { - if (logTerse()) { - log << "WARNING: All model parameters are fixed!"; - log << std::endl; - log << " ctlike will proceed without fitting parameters."; - log << std::endl; - log << " All curvature matrix elements will be zero."; - log << std::endl; - } + log_string(TERSE, "WARNING: All model parameters are fixed!"); + log_string(TERSE, " ctlike will proceed without fitting parameters."); + log_string(TERSE, " All curvature matrix elements will be zero."); } // Perform LM optimization @@ -494,11 +460,8 @@ void ctlike::optimize_lm(void) // Dump new header log.indent(0); - if (logTerse()) { - log << std::endl; - log.header1("Maximum likelihood re-optimisation"); - log.indent(1); - } + log_header1(TERSE, "Maximum likelihood re-optimisation"); + log.indent(1); // Optimise again m_obs.optimize(*m_opt); @@ -506,13 +469,10 @@ void ctlike::optimize_lm(void) } // Optionally show curvature matrix - if (logNormal()) { - log << std::endl; - log.header1("Curvature matrix"); - log.indent(1); - log << *(const_cast(m_obs.function()).curvature()); - log << std::endl; - } + log_header1(EXPLICIT, "Curvature matrix"); + log.indent(1); + log_string(EXPLICIT, (const_cast + (m_obs.function()).curvature())->print()); // Compute errors m_obs.errors(*m_opt); @@ -537,11 +497,8 @@ void ctlike::optimize_lm(void) double ctlike::reoptimize_lm(void) { // Write Header for optimization and indent for optimizer logging - if (logTerse()) { - log << std::endl; - log.header1("Maximum likelihood re-optimisation"); - log.indent(1); - } + log_header1(TERSE, "Maximum likelihood re-optimisation"); + log.indent(1); // Create a clone of the optimizer for the re-optimisation GOptimizer* opt = m_opt->clone(); @@ -559,11 +516,8 @@ double ctlike::reoptimize_lm(void) // Write optimization results log.indent(0); - if (logTerse()) { - log << std::endl; - log.header1("Maximum likelihood re-optimisation results"); - log << *opt << std::endl; - } + log_header1(NORMAL, "Maximum likelihood re-optimisation results"); + log_string(NORMAL, opt->print(m_chatter)); // Return return (logL); diff --git a/src/support/ctool.cpp b/src/support/ctool.cpp index 9eb8a3e..c66c6f1 100644 --- a/src/support/ctool.cpp +++ b/src/support/ctool.cpp @@ -52,6 +52,8 @@ #define G_SETUP_OBSERVATION "ctool::setup_observations(GObservations&)" #define G_GET_MEAN_POINTING "ctool::get_mean_pointing(GObservations&)" #define G_CREATE_EBOUNDS "ctool::create_ebounds()" +#define G_RESTORE_EDISP "ctool::restore_edisp(GObservations&, "\ + "std::vector&)" #define G_SET_OBS_RESPONSE "ctool::set_obs_response(GCTAObservation*)" #define G_PROVIDE_HELP "ctool::provide_help()" @@ -994,7 +996,7 @@ void ctool::set_response(GObservations& obs) // Is this observation a CTA observation? GCTAObservation* cta = dynamic_cast(obs[i]); - // Yes ... + // Yes, the set the response if the observation does not yet have one if (cta != NULL) { // Set response if we don't have one @@ -1012,6 +1014,87 @@ void ctool::set_response(GObservations& obs) /***********************************************************************//** + * @brief Set energy dispersion to CTA observations + * + * @param[in] obs Observation container. + * @param[in] edisp Requested energy dispersion flag. + * @return Vector of old energy dispersion flags. + * + * Applies energy dispersion to all CTA observations. The method returns a + * vector with the old energy dispersion flags. + ***************************************************************************/ +std::vector ctool::set_edisp(GObservations& obs, const bool& edisp) const +{ + // Allocate vector for energy dispersion flags of all observation in + // the container + std::vector old_edisp(obs.size(), false); + + // Loop over all observations in observation container + for (int i = 0; i < obs.size(); ++i) { + + // Is this observation a CTA observation? + GCTAObservation* cta = dynamic_cast(obs[i]); + + // Yes, then set the energy dispersion flag + if (cta != NULL) { + + // Save old energy dispersion flag + old_edisp[i] = cta->response()->apply_edisp(); + + // Set energy dispersion flag according to input parameter + cta->response()->apply_edisp(edisp); + + } // endif: observation was CTA observation + + } // endfor: loop over all observations in container + + // Return old energy dispersion flags + return old_edisp; +} + + +/***********************************************************************//** + * @brief Restore energy dispersion flags of CTA observations + * + * @param[in] obs Observation container. + * @param[in] edisp Vector of energy dispersion flags. + * + * Restores energy dispersion flags that have previously been returned by + * the set_edisp() method. The number of observations between the calls + * of the set_edisp() and this method has to be the same. + ***************************************************************************/ +void ctool::restore_edisp(GObservations& obs, const std::vector& edisp) const +{ + // Check that energy dispersion flag vector has the same size as the + // observation container + if (edisp.size() != obs.size()) { + std::string msg = "Size of energy dispersion flag vector ("+ + gammalib::str(edisp.size())+" elements) is " + "incompatible with number of observations in the " + "observation container ("+gammalib::str(obs.size())+ + " observation)."; + throw GException::invalid_value(G_RESTORE_EDISP, msg); + } + + // Loop over all observations in observation container + for (int i = 0; i < obs.size(); ++i) { + + // Is this observation a CTA observation? + GCTAObservation* cta = dynamic_cast(obs[i]); + + // Yes, then restore the energy dispersion flag + if (cta != NULL) { + cta->response()->apply_edisp(edisp[i]); + } + + } // endfor: loop over all observations in container + + // Return + return; +} + + +/***********************************************************************//** * @brief Set response for CTA observation * * @param[in,out] obs CTA observation diff --git a/src/support/ctool.hpp b/src/support/ctool.hpp index b2dabf0..7178f0a 100644 --- a/src/support/ctool.hpp +++ b/src/support/ctool.hpp @@ -94,14 +94,17 @@ protected: const std::string& what = "Model"); // Protected support methods - void set_response(GObservations& obs); - void set_obs_response(GCTAObservation* obs); - void set_obs_bounds(GObservations& obs); - GSkyDir get_mean_pointing(const GObservations& obs); - size_t get_current_rss(void); - std::string get_obs_header(const GObservation* obs); - GEnergies insert_energy_boundaries(const GEnergies& energies, - const GCTAObservation& obs); + void set_response(GObservations& obs); + std::vector set_edisp(GObservations& obs, const bool& edisp) const; + void restore_edisp(GObservations& obs, + const std::vector& edisp) const; + void set_obs_response(GCTAObservation* obs); + void set_obs_bounds(GObservations& obs); + GSkyDir get_mean_pointing(const GObservations& obs); + size_t get_current_rss(void); + std::string get_obs_header(const GObservation* obs); + GEnergies insert_energy_boundaries(const GEnergies& energies, + const GCTAObservation& obs); // Protected warning strings std::string warn_too_few_energies(const GEnergies& energies) const; diff --git a/test/test_cterror.py b/test/test_cterror.py index cfd563e..da62a1c 100755 --- a/test/test_cterror.py +++ b/test/test_cterror.py @@ -101,14 +101,31 @@ class Test(test): """ Test cterror from Python """ - # Set-up cterror + # Allocate cterror error = ctools.cterror() + + # Check that empty cterror tool holds an empty observation and + # optimizer + self._check_obs(error.obs(), nobs=0, nmodels=0) + self._check_opt(error.opt()) + + # Check that saving saves an empty model definition file + error['outmodel'] = 'cterror_py0.xml' + error['logfile'] = 'cterror_py0.log' + error.logFileOpen() + error.save() + self._check_result_file('cterror_py0.xml', nmodels=0) + + # Check that clearing does not lead to an exception or segfault + error.clear() + + # Now set cterror parameters error['inobs'] = self._events error['inmodel'] = self._model - error['outmodel'] = 'cterror_py1.xml' error['srcname'] = 'Crab' error['caldb'] = self._caldb error['irf'] = self._irf + error['outmodel'] = 'cterror_py1.xml' error['logfile'] = 'cterror_py1.log' error['chatter'] = 2 @@ -120,38 +137,130 @@ class Test(test): # Check result file self._check_result_file('cterror_py1.xml') + # Copy cterror tool + cpy_error = error.copy() + + # Check observation and optimizer of copy + self._check_obs(cpy_error.obs()) + self._check_opt(cpy_error.opt()) + + # Execute copy of cterror tool again, now with a higher chatter + # level than before + cpy_error['outmodel'] = 'cterror_py2.xml' + cpy_error['logfile'] = 'cterror_py2.log' + cpy_error['chatter'] = 3 + cpy_error.logFileOpen() # Needed to get a new log file + cpy_error.execute() + + # Check result file + self._check_result_file('cterror_py2.xml') + + # Now clear copy of cterror tool + cpy_error.clear() + + # Check that the cleared copy has also cleared the observations + self._check_obs(cpy_error.obs(), nobs=0, nmodels=0) + # Return return # Check result file - def _check_result_file(self, filename): + def _check_result_file(self, filename, nmodels=2): """ Check result file + + Parameters + ---------- + filename : str + Model definition XML file + nmodels : int, optional + Expected number of models + """ + # Open models from result file + models = gammalib.GModels(filename) + + # Check models + self._check_models(models, nmodels=nmodels) + + # Return + return + + # Check observation and models + def _check_obs(self, obs, nobs=1, nmodels=2): + """ + Check observation and models + + Parameters + ---------- + obs : `~gammalib.GObservations` + Models + nobs : int, optional + Expected number of observations + nmodels : int, optional + Expected number of models + """ + # Check number of observations + self.test_value(obs.size(), nobs, 'Check number of observations') + + # Check models + self._check_models(obs.models(), nmodels=nmodels) + + # Return + return + + # Check models + def _check_models(self, models, nmodels=2): + """ + Check observation and models + + Parameters + ---------- + models : `~gammalib.GModels` + Models + nmodels : int, optional + Expected number of models + """ + # Check number of models + self.test_value(models.size(), nmodels, 'Check number of models') + + # If there are models in the container then check them + if nmodels > 0: + self.test_value(models['Crab']['Prefactor'].value(), 1.58907e-16, + 1.0e-3, 'Check fitted Crab Prefactor') + self.test_value(models['Crab']['Prefactor'].error(), 0.0529105e-16, + 1.0e-3, 'Check Crab Prefactor error') + self.test_value(models['Crab']['Index'].value(), -2.43549, + 1.0e-3, 'Check fitted Crab Index') + self.test_value(models['Crab']['Index'].error(), 0.027804, + 1.0e-3, 'Check Crab Index error') + self.test_value(models['Background']['Prefactor'].value(), 61.6919e-6, + 1.0e-3, 'Check fitted background Prefactor') + self.test_value(models['Background']['Prefactor'].error(), 1.49438e-6, + 1.0e-3, 'Check background Prefactor error') + self.test_value(models['Background']['Index'].value(), -2.20535, + 1.0e-3, 'Check fitted background Index') + self.test_value(models['Background']['Index'].error(), 0.0113269, + 1.0e-3, 'Check background Index error') + self.test_value(models['Background']['Sigma'].value(), 3.04252, + 1.0e-3, 'Check fitted background Sigma') + self.test_value(models['Background']['Sigma'].error(), 0.0307008, + 1.0e-3, 'Check background Sigma error') + + # Return + return + + # Check optimizer + def _check_opt(self, opt): + """ + Check optimizer + + Parameters + ---------- + opt : `~gammalib.GOptimizerLM` + Optimizer """ - # Open result file - result = gammalib.GModels(filename) - - # Check results - self.test_value(result['Crab']['Prefactor'].value(), 1.58907e-16, - 1.0e-3, 'Check fitted Crab Prefactor') - self.test_value(result['Crab']['Prefactor'].error(), 0.0529105e-16, - 1.0e-3, 'Check Crab Prefactor error') - self.test_value(result['Crab']['Index'].value(), -2.43549, - 1.0e-3, 'Check fitted Crab Index') - self.test_value(result['Crab']['Index'].error(), 0.027804, - 1.0e-3, 'Check Crab Index error') - self.test_value(result['Background']['Prefactor'].value(), 61.6919e-6, - 1.0e-3, 'Check fitted background Prefactor') - self.test_value(result['Background']['Prefactor'].error(), 1.49438e-6, - 1.0e-3, 'Check background Prefactor error') - self.test_value(result['Background']['Index'].value(), -2.20535, - 1.0e-3, 'Check fitted background Index') - self.test_value(result['Background']['Index'].error(), 0.0113269, - 1.0e-3, 'Check background Index error') - self.test_value(result['Background']['Sigma'].value(), 3.04252, - 1.0e-3, 'Check fitted background Sigma') - self.test_value(result['Background']['Sigma'].error(), 0.0307008, - 1.0e-3, 'Check background Sigma error') + # Check optimizer + self.test_value(opt.status(), 0, 'Check optimizer status') # Return return diff --git a/test/test_ctlike.py b/test/test_ctlike.py index 77c3e68..932f853 100755 --- a/test/test_ctlike.py +++ b/test/test_ctlike.py @@ -101,8 +101,25 @@ class Test(test): """ Test ctlike from Python """ - # Set-up ctlike + # Allocate ctlike like = ctools.ctlike() + + # Check that empty ctlike tool holds an empty observation and + # optimizer + self._check_obs(like.obs(), nobs=0, nmodels=0) + self._check_opt(like.opt()) + + # Check that saving saves an empty model definition file + like['outmodel'] = 'ctlike_py0.xml' + like['logfile'] = 'ctlike_py0.log' + like.logFileOpen() + like.save() + self._check_result_file('ctlike_py0.xml', nmodels=0) + + # Check that clearing does not lead to an exception or segfault + like.clear() + + # Now set ctlike parameters like['inobs'] = self._events like['inmodel'] = self._model like['caldb'] = self._caldb @@ -119,46 +136,130 @@ class Test(test): # Check result file self._check_result_file('ctlike_py1.xml') - # Retrieve observation container, save and check model file - obs = like.obs() - obs.models().save('ctlike_py2.xml') + # Copy ctlike tool + cpy_like = like.copy() + + # Check observation and optimizer of copy + self._check_obs(cpy_like.obs()) + self._check_opt(cpy_like.opt()) + + # Execute copy of ctlike tool again, now with a higher chatter + # level than before + cpy_like['outmodel'] = 'ctlike_py2.xml' + cpy_like['logfile'] = 'ctlike_py2.log' + cpy_like['chatter'] = 3 + cpy_like.logFileOpen() # Needed to get a new log file + cpy_like.execute() + + # Check result file self._check_result_file('ctlike_py2.xml') - # Retrieve optimizer - opt = like.opt() + # Now clear copy of ctlike tool + cpy_like.clear() + + # Check that the cleared copy has also cleared the observations + self._check_obs(cpy_like.obs(), nobs=0, nmodels=0) # Return return # Check result file - def _check_result_file(self, filename): + def _check_result_file(self, filename, nmodels=2): """ Check result file + + Parameters + ---------- + filename : str + Model definition XML file + nmodels : int, optional + Expected number of models + """ + # Open models from result file + models = gammalib.GModels(filename) + + # Check models + self._check_models(models, nmodels=nmodels) + + # Return + return + + # Check observation and models + def _check_obs(self, obs, nobs=1, nmodels=2): + """ + Check observation and models + + Parameters + ---------- + obs : `~gammalib.GObservations` + Models + nobs : int, optional + Expected number of observations + nmodels : int, optional + Expected number of models + """ + # Check number of observations + self.test_value(obs.size(), nobs, 'Check number of observations') + + # Check models + self._check_models(obs.models(), nmodels=nmodels) + + # Return + return + + # Check models + def _check_models(self, models, nmodels=2): + """ + Check observation and models + + Parameters + ---------- + models : `~gammalib.GModels` + Models + nmodels : int, optional + Expected number of models + """ + # Check number of models + self.test_value(models.size(), nmodels, 'Check number of models') + + # If there are models in the container then check them + if nmodels > 0: + self.test_value(models['Crab']['Prefactor'].value(), 1.58907e-16, + 1.0e-3, 'Check fitted Crab Prefactor') + self.test_value(models['Crab']['Prefactor'].error(), 0.0526982e-16, + 1.0e-3, 'Check Crab Prefactor error') + self.test_value(models['Crab']['Index'].value(), -2.43549, + 1.0e-3, 'Check fitted Crab Index') + self.test_value(models['Crab']['Index'].error(), 0.0248116, + 1.0e-3, 'Check Crab Index error') + self.test_value(models['Background']['Prefactor'].value(), 61.6919e-6, + 1.0e-3, 'Check fitted background Prefactor') + self.test_value(models['Background']['Prefactor'].error(), 1.49438e-6, + 1.0e-3, 'Check background Prefactor error') + self.test_value(models['Background']['Index'].value(), -2.20535, + 1.0e-3, 'Check fitted background Index') + self.test_value(models['Background']['Index'].error(), 0.0113269, + 1.0e-3, 'Check background Index error') + self.test_value(models['Background']['Sigma'].value(), 3.04252, + 1.0e-3, 'Check fitted background Sigma') + self.test_value(models['Background']['Sigma'].error(), 0.0307008, + 1.0e-3, 'Check background Sigma error') + + # Return + return + + # Check optimizer + def _check_opt(self, opt): + """ + Check optimizer + + Parameters + ---------- + opt : `~gammalib.GOptimizerLM` + Optimizer """ - # Open result file - result = gammalib.GModels(filename) - - # Check results - self.test_value(result['Crab']['Prefactor'].value(), 1.58907e-16, - 1.0e-3, 'Check fitted Crab Prefactor') - self.test_value(result['Crab']['Prefactor'].error(), 0.0526982e-16, - 1.0e-3, 'Check Crab Prefactor error') - self.test_value(result['Crab']['Index'].value(), -2.43549, - 1.0e-3, 'Check fitted Crab Index') - self.test_value(result['Crab']['Index'].error(), 0.0248116, - 1.0e-3, 'Check Crab Index error') - self.test_value(result['Background']['Prefactor'].value(), 61.6919e-6, - 1.0e-3, 'Check fitted background Prefactor') - self.test_value(result['Background']['Prefactor'].error(), 1.49438e-6, - 1.0e-3, 'Check background Prefactor error') - self.test_value(result['Background']['Index'].value(), -2.20535, - 1.0e-3, 'Check fitted background Index') - self.test_value(result['Background']['Index'].error(), 0.0113269, - 1.0e-3, 'Check background Index error') - self.test_value(result['Background']['Sigma'].value(), 3.04252, - 1.0e-3, 'Check fitted background Sigma') - self.test_value(result['Background']['Sigma'].error(), 0.0307008, - 1.0e-3, 'Check background Sigma error') + # Check optimizer + self.test_value(opt.status(), 0, 'Check optimizer status') # Return return