From 13d498b918b48ba868c1e3d680f78bd4295138a2 Mon Sep 17 00:00:00 2001 From: knodlseder <> Date: Fri, 12 Nov 2010 21:54:03 +0000 Subject: [PATCH] generalize code and implement logger interface --- src/ctlike/ctlike.cpp | 203 ++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++------------------------------------------------------------- src/ctlike/ctlike.hpp | 17 +++++++++++++---- 2 files changed, 155 insertions(+), 65 deletions(-) diff --git a/src/ctlike/ctlike.cpp b/src/ctlike/ctlike.cpp index 8dab2b0..3069337 100644 --- a/src/ctlike/ctlike.cpp +++ b/src/ctlike/ctlike.cpp @@ -34,11 +34,14 @@ /***********************************************************************//** * @brief Void constructor ***************************************************************************/ -ctlike::ctlike(void) : GApplication(CTLIKE_NAME, CTLIKE_VERSION) +ctlike::ctlike(void) : GApplication(CTLIKE_NAME, CTLIKE_VERSION) { // Initialise members init_members(); + // Write header into logger + log_header(); + // Return return; } @@ -56,6 +59,9 @@ ctlike::ctlike(int argc, char *argv[]) : // Initialise members init_members(); + // Write header into logger + log_header(); + // Return return; } @@ -88,75 +94,165 @@ ctlike::~ctlike(void) /***********************************************************************//** * @brief Run maximum likelihood analysis + * + * The following analysis steps are performed: + * 1. Read the parameters (and write them into logger) + * 2. Load observation + * 3. Setup models for optimizing + * 4. Optimize model (and write result into logger) + * 5. Write results out as XML model ***************************************************************************/ -int ctlike::run(void) +void ctlike::run(void) { - // Initialise return code - int rc = 0; + // Get parameters + get_parameters(); + + // Write parameters into logger + if (logTerse()) { + log_parameters(); + log << std::endl; + } + + // Load observation + if (m_method == "UNBINNED") + unbinned(m_evfile); + else + binned(m_cntmap); + + // Setup models for optimizing. + m_models = GModels(m_srcmdl); + m_obs.models(m_models); + + // Optimize + optimize_lm(); + + // Write results into logger + if (logTerse()) { + log << " Maximum log likelihood ....: " << m_logL << std::endl; + log << " Npred .....................: " << m_obs.npred() << std::endl; + log << m_models << std::endl << std::endl; + } + + // Write results out as XML model + if (toupper(m_outmdl) != "NONE") + m_models.save(m_outmdl); + + // Return + return; +} + +/***********************************************************************//** + * @brief Get application parameters + * + * Get all task parameters from parameter file or (if required) by querying + * the user. + ***************************************************************************/ +void ctlike::get_parameters(void) +{ // Get standard parameters m_method = toupper(par("method")->value()); m_caldb = par("caldb")->value(); m_irf = par("irf")->value(); m_srcmdl = par("srcmdl")->value(); m_outmdl = par("outmdl")->value(); - m_models = GModels(m_srcmdl); - // Set fixed parameters - m_opt = new GOptimizerLM; + // Get unbinned parameters ... + if (m_method == "UNBINNED") { + m_evfile = par("evfile")->value(); + } + + // ... or get binned parameters + else if (m_method == "BINNED") { + m_cntmap = par("cntmap")->value(); + } + + // Return + return; +} + + +/***********************************************************************//** + * @brief Optimize using Levenberg-Marquardt method + ***************************************************************************/ +void ctlike::optimize_lm(void) +{ + // Allocate optimizer. The logger is only passed to the optimizer + // constructor if optimizer logging is requested. + GOptimizerLM* opt = (logExplicit()) ? new GOptimizerLM(log) + : new GOptimizerLM(); + + // Assign optimizer + m_opt = opt; // Set optimizer parameters - ((GOptimizerLM*)m_opt)->max_iter(m_max_iter); + opt->max_iter(m_max_iter); - // Branch on method - if (m_method == "UNBINNED") - rc = unbinned(); - else - rc = binned(); + // Write Header for optimization and indent for optimizer logging + if (logExplicit()) { + log.header1("Maximum likelihood optimisation"); + log.indent(1); + } + + // Perform LM optimization + m_obs.optimize(*m_opt); - // Write result as XML model - if (toupper(m_outmdl) != "NONE") - m_models.save(m_outmdl); + // Get models back + m_models = *(m_obs.models()); + + // Store maximum log likelihood value + m_logL = -opt->value(); - // Dump results - std::cout << m_models << std::endl; + // Write optimization results + log.indent(0); + if (logTerse()) { + log << std::endl; + log.header1("Maximum likelihood optimisation results"); + log << *opt << std::endl; + } // Free optimizer - delete m_opt; + delete opt; + + // Reset optimizer pointer + m_opt = NULL; // Return - return rc; + return; } /***********************************************************************//** - * @brief Perform unbinned maximum likelihood analysis + * @brief Load unbinned observation + * + * @param[in] evfile Events file name. + * + * @todo ROI, energy range and time range information is hardwired. This is + * not needed in the future once this information is correctly extracted + * from the event file header. ***************************************************************************/ -int ctlike::unbinned(void) +void ctlike::unbinned(const std::string& evfile) { - // Initialise return code - int rc = 0; - - // Get parameters - std::string evfile = par("evfile")->value(); - // Declare observations GCTAObservation run; - // Setup ROI covered by data + // DUMMY: Setup ROI covered by data. This is needed since actual test + // data do not contain any ROI informations. GCTAInstDir instDir; GCTARoi roi; instDir.radec_deg(117.02, -33.35); // Adapt to file roi.centre(instDir); roi.radius(2.5); - // Setup energy range covered by data + // DUMMY: Setup energy range covered by data. This is needed since actual + // test data do not contain any energy range information GEnergy emin; GEnergy emax; emin.TeV(0.02); emax.TeV(100.0); - // Setup time range covered by data + // DUMMY: Setup time range covered by data. This is needed since actual + // test data do not contain any time range information GTime tstart; GTime tstop; tstart.met(0.0); @@ -171,50 +267,32 @@ int ctlike::unbinned(void) run.ebounds()->append(emin, emax); m_obs.append(run); - // Setup models for optimizing - m_obs.models(m_models); - - // Perform LM optimization - m_obs.optimize(*m_opt); - - // Get models back - m_models = *(m_obs.models()); - // Return - return rc; + return; } /***********************************************************************//** - * @brief Perform binned maximum likelihood analysis + * @brief Load binned observation + * + * @param[in] cntmap Counts map file name. ***************************************************************************/ -int ctlike::binned(void) +void ctlike::binned(const std::string& cntmap) { - // Initialise return code - int rc = 0; - - // Get parameters - std::string cntmap = par("cntmap")->value(); - // Declare observations GCTAObservation run; // Load binned CTA observation run.load_binned(cntmap); - run.response(m_irf, m_caldb); - m_obs.append(run); - // Setup models for optimizing - m_obs.models(m_models); - - // Perform LM optimization - m_obs.optimize(*m_opt); + // Set response + run.response(m_irf, m_caldb); - // Get models back - m_models = *(m_obs.models()); + // Append observation to container + m_obs.append(run); // Return - return rc; + return; } @@ -236,10 +314,14 @@ void ctlike::init_members(void) m_srcmdl.clear(); m_outmdl.clear(); m_models.clear(); - //m_obs.clear(); //!< NOT YET IMPLEMENTED + m_obs.clear(); m_max_iter = 1000; + m_logL = 0.0; m_opt = NULL; + // Set logger properties + log.date(true); + // Return return; } @@ -251,7 +333,6 @@ void ctlike::init_members(void) void ctlike::free_members(void) { // Free members - delete m_opt; // Return return; diff --git a/src/ctlike/ctlike.hpp b/src/ctlike/ctlike.hpp index 2cc7ab0..32008be 100644 --- a/src/ctlike/ctlike.hpp +++ b/src/ctlike/ctlike.hpp @@ -42,9 +42,11 @@ public: ~ctlike(void); // Methods - int run(void); - int unbinned(void); - int binned(void); + void run(void); + void get_parameters(void); + void optimize_lm(void); + void unbinned(const std::string& evfile); + void binned(const std::string& cntmap); protected: // Protected methods @@ -60,8 +62,15 @@ protected: GModels m_models; //!< Source models GObservations m_obs; //!< Observations - // Fixed parameters + // Unbinned likelihood parameters + std::string m_evfile; //!< Events file + + // Binned likelihood parameters + std::string m_cntmap; //!< Counts map + + // Internal parameters int m_max_iter; //!< Maximum number of iterations + double m_logL; //!< Maximum log likelihood GOptimizer* m_opt; //!< Optimizer }; -- libgit2 0.21.4