Lili Yang / ctools

Commit e9ae0dd9457fe4d86a8665c83799b3357fa38334

Authored by knodlseder
1 parent 5e271e1d

update to support standardized Python interface

Inline Side-by-side
Showing 4 changed files with 309 additions and 89 deletions
pyext/ctlike.i
  View file @e9ae0dd
1 1 /***************************************************************************
2 2 * ctlike - CTA maximum likelihood tool *
3 3 * ----------------------------------------------------------------------- *
4   - * copyright (C) 2010 by Jurgen Knodlseder *
  4 + * copyright (C) 2010-2011 by Jurgen Knodlseder *
5 5 * ----------------------------------------------------------------------- *
6 6 * *
7 7 * This program is free software; you can redistribute it and/or modify *
... ... @@ -12,7 +12,7 @@
12 12 ***************************************************************************/
13 13 /**
14 14 * @file ctlike.i
15   - * @brief CTA maximum likelihood tool SWIG definition
  15 + * @brief CTA maximum likelihood tool Python interface definition
16 16 * @author J. Knodlseder
17 17 */
18 18 %{
... ... @@ -25,21 +25,27 @@
25 25 /***********************************************************************//**
26 26 * @class ctlike
27 27 *
28   - * @brief CTA maximum likelihood tool SWIG interface defintion.
  28 + * @brief CTA maximum likelihood tool Python interface defintion
29 29 ***************************************************************************/
30 30 class ctlike : public GApplication {
31 31 public:
32 32 // Constructors and destructors
33 33 ctlike(void);
  34 + explicit ctlike(GObservations obs);
34 35 ctlike(int argc, char *argv[]);
35   - ~ctlike(void);
  36 + ctlike(const ctlike& app);
  37 + virtual ~ctlike(void);
36 38  
37 39 // Methods
38   - void run(void);
39   - void get_parameters(void);
40   - void optimize_lm(void);
41   - void unbinned(const std::string& evfile);
42   - void binned(const std::string& cntmap);
  40 + void clear(void);
  41 + void execute(void);
  42 + void run(void);
  43 + void save(void);
  44 + GObservations& obs(void) { return m_obs; }
  45 + void get_parameters(void);
  46 + void optimize_lm(void);
  47 + void load_unbinned(const std::string& evfile);
  48 + void load_binned(const std::string& cntmap);
43 49 };
44 50  
45 51  
... ...
src/ctlike/ctlike.cpp
  View file @e9ae0dd
... ... @@ -24,6 +24,13 @@
24 24 #include "ctlike.hpp"
25 25 #include "GTools.hpp"
26 26  
  27 +/* __ Method name definitions ____________________________________________ */
  28 +#define G_GET_PARAMETERS "ctlike::get_parameters()"
  29 +
  30 +/* __ Debug definitions __________________________________________________ */
  31 +
  32 +/* __ Coding definitions _________________________________________________ */
  33 +
27 34  
28 35 /*==========================================================================
29 36 = =
... ... @@ -48,6 +55,28 @@ ctlike::ctlike(void) : GApplication(CTLIKE_NAME, CTLIKE_VERSION)
48 55  
49 56  
50 57 /***********************************************************************//**
  58 + * @brief Observations constructor
  59 + *
  60 + * This method creates an instance of the class by copying an existing
  61 + * observations container.
  62 + ***************************************************************************/
  63 +ctlike::ctlike(GObservations obs) : GApplication(CTLIKE_NAME, CTLIKE_VERSION)
  64 +{
  65 + // Initialise members
  66 + init_members();
  67 +
  68 + // Set observations
  69 + m_obs = obs;
  70 +
  71 + // Write header into logger
  72 + log_header();
  73 +
  74 + // Return
  75 + return;
  76 +}
  77 +
  78 +
  79 +/***********************************************************************//**
51 80 * @brief Command line constructor
52 81 *
53 82 * @param[in] argc Number of arguments in command line.
... ... @@ -68,6 +97,24 @@ ctlike::ctlike(int argc, char *argv[]) :
68 97  
69 98  
70 99 /***********************************************************************//**
  100 + * @brief Copy constructor
  101 + *
  102 + * @param[in] app Application.
  103 + ***************************************************************************/
  104 +ctlike::ctlike(const ctlike& app) : GApplication(app)
  105 +{
  106 + // Initialise members
  107 + init_members();
  108 +
  109 + // Copy members
  110 + copy_members(app);
  111 +
  112 + // Return
  113 + return;
  114 +}
  115 +
  116 +
  117 +/***********************************************************************//**
71 118 * @brief Destructor
72 119 ***************************************************************************/
73 120 ctlike::~ctlike(void)
... ... @@ -86,6 +133,35 @@ ctlike::~ctlike(void)
86 133 = =
87 134 ==========================================================================*/
88 135  
  136 +/***********************************************************************//**
  137 + * @brief Assignment operator
  138 + *
  139 + * @param[in] app Application.
  140 + ***************************************************************************/
  141 +ctlike& ctlike::operator= (const ctlike& app)
  142 +{
  143 + // Execute only if object is not identical
  144 + if (this != &app) {
  145 +
  146 + // Copy base class members
  147 + this->GApplication::operator=(app);
  148 +
  149 + // Free members
  150 + free_members();
  151 +
  152 + // Initialise members
  153 + init_members();
  154 +
  155 + // Copy members
  156 + copy_members(app);
  157 +
  158 + } // endif: object was not identical
  159 +
  160 + // Return this object
  161 + return *this;
  162 +}
  163 +
  164 +
89 165 /*==========================================================================
90 166 = =
91 167 = Public methods =
... ... @@ -93,6 +169,47 @@ ctlike::~ctlike(void)
93 169 ==========================================================================*/
94 170  
95 171 /***********************************************************************//**
  172 + * @brief Clear instance
  173 + ***************************************************************************/
  174 +void ctlike::clear(void)
  175 +{
  176 + // Free members
  177 + free_members();
  178 + this->GApplication::free_members();
  179 +
  180 + // Initialise members
  181 + this->GApplication::init_members();
  182 + init_members();
  183 +
  184 + // Return
  185 + return;
  186 +}
  187 +
  188 +
  189 +/***********************************************************************//**
  190 + * @brief Execute application
  191 + *
  192 + * This method performs a maximum likelihood analysis of a observation given
  193 + * in an observation container and saves the results in an XML file.
  194 + ***************************************************************************/
  195 +void ctlike::execute(void)
  196 +{
  197 + // Read ahead output filename so that it gets dumped correctly in the
  198 + // parameters log
  199 + m_outmdl = par("outmdl")->value();
  200 +
  201 + // Bin the event data
  202 + run();
  203 +
  204 + // Save the results into XML file
  205 + save();
  206 +
  207 + // Return
  208 + return;
  209 +}
  210 +
  211 +
  212 +/***********************************************************************//**
96 213 * @brief Run maximum likelihood analysis
97 214 *
98 215 * The following analysis steps are performed:
... ... @@ -100,7 +217,6 @@ ctlike::~ctlike(void)
100 217 * 2. Load observation
101 218 * 3. Setup models for optimizing
102 219 * 4. Optimize model (and write result into logger)
103   - * 5. Write results out as XML model
104 220 ***************************************************************************/
105 221 void ctlike::run(void)
106 222 {
... ... @@ -117,36 +233,50 @@ void ctlike::run(void)
117 233 log << std::endl;
118 234 }
119 235  
120   - // Load observation
121   - if (m_method == "UNBINNED")
122   - unbinned(m_evfile);
123   - else
124   - binned(m_cntmap);
125   -
126   - // Write observation into logger
  236 + // Write observation(s) into logger
127 237 if (logTerse()) {
128 238 log << std::endl;
129   - log.header1("Observation");
  239 + if (m_obs.size() > 1)
  240 + log.header1("Observations");
  241 + else
  242 + log.header1("Observation");
130 243 log << m_obs << std::endl;
131 244 }
132 245  
133   - // Setup models for optimizing.
134   - m_models = GModels(m_srcmdl);
135   - m_obs.models(m_models);
136   -
137   - // Optimize
  246 + // Optimize model parameters using LM optimizer
138 247 optimize_lm();
139 248  
140 249 // Write results into logger
141 250 if (logTerse()) {
142 251 log << " Maximum log likelihood ....: " << m_logL << std::endl;
143 252 log << " Npred .....................: " << m_obs.npred() << std::endl;
144   - log << m_models << std::endl << std::endl;
  253 + log << m_obs.models() << std::endl;
145 254 }
146   -
  255 +
  256 + // Return
  257 + return;
  258 +}
  259 +
  260 +
  261 +/***********************************************************************//**
  262 + * @brief Save results
  263 + *
  264 + * This method saves the fit results in a XML file.
  265 + ***************************************************************************/
  266 +void ctlike::save(void)
  267 +{
  268 + // Write header
  269 + if (logTerse()) {
  270 + log << std::endl;
  271 + log.header1("Save results");
  272 + }
  273 +
  274 + // Get output filename
  275 + m_outmdl = par("outmdl")->value();
  276 +
147 277 // Write results out as XML model
148 278 if (toupper(m_outmdl) != "NONE")
149   - m_models.save(m_outmdl);
  279 + m_obs.models().save(m_outmdl);
150 280  
151 281 // Return
152 282 return;
... ... @@ -156,29 +286,78 @@ void ctlike::run(void)
156 286 /***********************************************************************//**
157 287 * @brief Get application parameters
158 288 *
159   - * Get all task parameters from parameter file or (if required) by querying
160   - * the user.
  289 + * GException::par_error
  290 + * Invalid analysis method specified (neither BINNED nor UNBINNED).
  291 + *
  292 + * Get all required task parameters from the parameter file or (if specified)
  293 + * by querying the user. Observation dependent parameters will only be read
  294 + * if the observation container is actually empty. Observation dependent
  295 + * parameters are:
  296 + * "method" (analysis method),
  297 + * "stat" (statistics to be used for observation),
  298 + * "caldb" (calibration database),
  299 + * "irf" (instrument response function),
  300 + * "evfile" (event file name), and
  301 + * "cntmap" (counts map file name).
  302 + * The model will only be loaded if no model components exist in the
  303 + * observation container.
161 304 ***************************************************************************/
162 305 void ctlike::get_parameters(void)
163 306 {
  307 + // If there is are no models associated with the observations then
  308 + // load now the model definition
  309 + if (m_obs.models().size() == 0) {
  310 +
  311 + // Get models XML filename
  312 + std::string filename = par("srcmdl")->value();
  313 +
  314 + // Setup models for optimizing.
  315 + m_obs.models(GModels(filename));
  316 +
  317 + } // endif: no models were associated with observations
  318 +
  319 + // If there are no observations in container then add a single CTA
  320 + // observation using the parameters from the parameter file
  321 + if (m_obs.size() == 0) {
  322 +
  323 + // Get observation parameters
  324 + m_method = toupper(par("method")->value());
  325 + m_stat = toupper(par("stat")->value());
  326 + m_caldb = par("caldb")->value();
  327 + m_irf = par("irf")->value();
  328 +
  329 + // Case A: set-up unbinned CTA observation
  330 + if (m_method == "UNBINNED") {
  331 +
  332 + // Get event file name
  333 + std::string filename = par("evfile")->value();
  334 +
  335 + // Load and append unbinned CTA observation
  336 + load_unbinned(filename);
  337 +
  338 + } // endif: unbinned analysis mode
  339 +
  340 + // Case B: set-up binned CTA observation
  341 + else if (m_method == "BINNED") {
  342 +
  343 + // Get counts map file name
  344 + std::string filename = par("cntmap")->value();
  345 +
  346 + // Load and append binned CTA observation
  347 + load_binned(filename);
  348 +
  349 + } // endif: binned analysis mode
  350 +
  351 + // ... otherwise signal an invalid analysis method
  352 + else {
  353 + throw GException::par_error(G_GET_PARAMETERS, "method",
  354 + "only \"BINNED\" or \"UNBINNED\" supported.");
  355 + }
  356 +
  357 + } // endif: there was no observation in the container
  358 +
164 359 // Get standard parameters
165   - m_method = toupper(par("method")->value());
166   - m_stat = toupper(par("stat")->value());
167 360 m_refit = par("refit")->boolean();
168   - m_caldb = par("caldb")->value();
169   - m_irf = par("irf")->value();
170   - m_srcmdl = par("srcmdl")->value();
171   - m_outmdl = par("outmdl")->value();
172   -
173   - // Get unbinned parameters ...
174   - if (m_method == "UNBINNED") {
175   - m_evfile = par("evfile")->value();
176   - }
177   -
178   - // ... or get binned parameters
179   - else if (m_method == "BINNED") {
180   - m_cntmap = par("cntmap")->value();
181   - }
182 361  
183 362 // Return
184 363 return;
... ... @@ -186,15 +365,19 @@ void ctlike::get_parameters(void)
186 365  
187 366  
188 367 /***********************************************************************//**
189   - * @brief Optimize using Levenberg-Marquardt method
  368 + * @brief Optimize model parameters using Levenberg-Marquardt method
190 369 ***************************************************************************/
191 370 void ctlike::optimize_lm(void)
192 371 {
  372 + // Free any existing optimizer
  373 + if (m_opt != NULL) delete m_opt;
  374 + m_opt = NULL;
  375 +
193 376 // Allocate optimizer. The logger is only passed to the optimizer
194 377 // constructor if optimizer logging is requested.
195   - GOptimizerLM* opt = (logTerse()) ? new GOptimizerLM(log)
  378 + GOptimizerLM* opt = (logTerse()) ? new GOptimizerLM(log)
196 379 : new GOptimizerLM();
197   -
  380 +
198 381 // Assign optimizer
199 382 m_opt = opt;
200 383  
... ... @@ -207,45 +390,36 @@ void ctlike::optimize_lm(void)
207 390 log.header1("Maximum likelihood optimisation");
208 391 log.indent(1);
209 392 }
210   -
  393 +
211 394 // Perform LM optimization
212   - m_obs.optimize(*m_opt);
  395 + m_obs.optimize(*opt);
213 396  
214 397 // Optionally refit
215 398 if (m_refit)
216   - m_obs.optimize(*m_opt);
217   -
218   - // Get models back
219   - m_models = m_obs.models();
  399 + m_obs.optimize(*opt);
220 400  
221 401 // Store maximum log likelihood value
222   - m_logL = -opt->value();
  402 + m_logL = -(opt->value());
223 403  
224 404 // Write optimization results
225 405 log.indent(0);
226 406 if (logTerse()) {
227 407 log << std::endl;
228   - log.header1("Maximum likelihood optimisation results");
  408 + log.header1("Maximum likelihood optimization results");
229 409 log << *opt << std::endl;
230 410 }
231 411  
232   - // Free optimizer
233   - delete opt;
234   -
235   - // Reset optimizer pointer
236   - m_opt = NULL;
237   -
238 412 // Return
239 413 return;
240 414 }
241 415  
242 416  
243 417 /***********************************************************************//**
244   - * @brief Load unbinned observation and append to observations contained
  418 + * @brief Load unbinned observation and append it to observations container
245 419 *
246 420 * @param[in] evfile Events file name.
247 421 ***************************************************************************/
248   -void ctlike::unbinned(const std::string& evfile)
  422 +void ctlike::load_unbinned(const std::string& evfile)
249 423 {
250 424 // Declare CTA observation
251 425 GCTAObservation obs;
... ... @@ -268,11 +442,11 @@ void ctlike::unbinned(const std::string&amp; evfile)
268 442  
269 443  
270 444 /***********************************************************************//**
271   - * @brief Load binned observation and append to observations contained
  445 + * @brief Load binned observation and append it to observations container
272 446 *
273 447 * @param[in] cntmap Counts map file name.
274 448 ***************************************************************************/
275   -void ctlike::binned(const std::string& cntmap)
  449 +void ctlike::load_binned(const std::string& cntmap)
276 450 {
277 451 // Declare CTA observation
278 452 GCTAObservation obs;
... ... @@ -307,12 +481,12 @@ void ctlike::init_members(void)
307 481 {
308 482 // Initialise members
309 483 m_method.clear();
  484 + m_stat.clear();
310 485 m_caldb.clear();
311 486 m_irf.clear();
312   - m_srcmdl.clear();
313 487 m_outmdl.clear();
314   - m_models.clear();
315 488 m_obs.clear();
  489 + m_refit = false;
316 490 m_max_iter = 100; // Set maximum number of iterations to 100
317 491 m_logL = 0.0;
318 492 m_opt = NULL;
... ... @@ -326,11 +500,43 @@ void ctlike::init_members(void)
326 500  
327 501  
328 502 /***********************************************************************//**
  503 + * @brief Copy class members
  504 + *
  505 + * @param[in] app Application.
  506 + ***************************************************************************/
  507 +void ctlike::copy_members(const ctlike& app)
  508 +{
  509 + // Copy attributes
  510 + m_method = app.m_method;
  511 + m_stat = app.m_stat;
  512 + m_refit = app.m_refit;
  513 + m_caldb = app.m_caldb;
  514 + m_irf = app.m_irf;
  515 + m_outmdl = app.m_outmdl;
  516 + m_obs = app.m_obs;
  517 + m_max_iter = app.m_max_iter;
  518 + m_logL = app.m_logL;
  519 + m_opt = app.m_opt->clone();
  520 +
  521 + // Return
  522 + return;
  523 +}
  524 +
  525 +
  526 +/***********************************************************************//**
329 527 * @brief Delete class members
330 528 ***************************************************************************/
331 529 void ctlike::free_members(void)
332 530 {
333 531 // Free members
  532 + if (m_opt != NULL) delete m_opt;
  533 +
  534 + // Mark pointers as free
  535 + m_opt = NULL;
  536 +
  537 + // Write separator into logger
  538 + if (logTerse())
  539 + log << std::endl;
334 540  
335 541 // Return
336 542 return;
... ...
src/ctlike/ctlike.hpp
  View file @e9ae0dd
... ... @@ -12,7 +12,7 @@
12 12 ***************************************************************************/
13 13 /**
14 14 * @file ctlike.hpp
15   - * @brief CTA maximum likelihood tool definition
  15 + * @brief CTA maximum likelihood tool interface definition
16 16 * @author J. Knodlseder
17 17 */
18 18  
... ... @@ -26,7 +26,7 @@
26 26  
27 27 /* __Definitions _________________________________________________________ */
28 28 #define CTLIKE_NAME "ctlike"
29   -#define CTLIKE_VERSION "00-01-00"
  29 +#define CTLIKE_VERSION "00-02-00"
30 30  
31 31  
32 32 /***********************************************************************//**
... ... @@ -38,19 +38,29 @@ class ctlike : public GApplication {
38 38 public:
39 39 // Constructors and destructors
40 40 ctlike(void);
  41 + explicit ctlike(GObservations obs);
41 42 ctlike(int argc, char *argv[]);
42   - ~ctlike(void);
  43 + ctlike(const ctlike& app);
  44 + virtual ~ctlike(void);
  45 +
  46 + // Operators
  47 + ctlike& operator= (const ctlike& app);
43 48  
44 49 // Methods
45   - void run(void);
46   - void get_parameters(void);
47   - void optimize_lm(void);
48   - void unbinned(const std::string& evfile);
49   - void binned(const std::string& cntmap);
  50 + void clear(void);
  51 + void execute(void);
  52 + void run(void);
  53 + void save(void);
  54 + GObservations& obs(void) { return m_obs; }
  55 + void get_parameters(void);
  56 + void optimize_lm(void);
  57 + void load_unbinned(const std::string& evfile);
  58 + void load_binned(const std::string& cntmap);
50 59  
51 60 protected:
52 61 // Protected methods
53 62 void init_members(void);
  63 + void copy_members(const ctlike& app);
54 64 void free_members(void);
55 65  
56 66 // User parameters
... ... @@ -59,18 +69,10 @@ protected:
59 69 bool m_refit; //!< Refitting
60 70 std::string m_caldb; //!< Calibration database
61 71 std::string m_irf; //!< Instrument response functions
62   - std::string m_srcmdl; //!< Source model XML file
63 72 std::string m_outmdl; //!< Source model output XML file
64   - GModels m_models; //!< Source models
65   - GObservations m_obs; //!< Observations
66   -
67   - // Unbinned likelihood parameters
68   - std::string m_evfile; //!< Events file
69 73  
70   - // Binned likelihood parameters
71   - std::string m_cntmap; //!< Counts map
72   -
73   - // Internal parameters
  74 + // Members
  75 + GObservations m_obs; //!< Observations
74 76 int m_max_iter; //!< Maximum number of iterations
75 77 double m_logL; //!< Maximum log likelihood
76 78 GOptimizer* m_opt; //!< Optimizer
... ...
src/ctlike/main.cpp
  View file @e9ae0dd
1 1 /***************************************************************************
2 2 * ctlike - CTA maximum likelihood tool main code *
3 3 * ----------------------------------------------------------------------- *
4   - * copyright (C) 2010 by Jurgen Knodlseder *
  4 + * copyright (C) 2010-2011 by Jurgen Knodlseder *
5 5 * ----------------------------------------------------------------------- *
6 6 * *
7 7 * This program is free software; you can redistribute it and/or modify *
... ... @@ -20,7 +20,7 @@
20 20 #ifdef HAVE_CONFIG_H
21 21 #include <config.h>
22 22 #endif
23   -#include <stdio.h>
  23 +#include <cstdio>
24 24 #include "ctlike.hpp"
25 25  
26 26  
... ... @@ -29,6 +29,9 @@
29 29 *
30 30 * @param[in] argc Number of arguments
31 31 * @param[in] argv Arguments
  32 + *
  33 + * This is the main entry point of the ctlike application. It allocates a
  34 + * ctlike object and executes the application.
32 35 ***************************************************************************/
33 36 int main (int argc, char *argv[])
34 37 {
... ... @@ -40,7 +43,10 @@ int main (int argc, char *argv[])
40 43  
41 44 // Run application
42 45 try {
43   - application.run();
  46 + // Execute application
  47 + application.execute();
  48 +
  49 + // Signal success
44 50 rc = 0;
45 51 }
46 52 catch (std::exception &e) {
... ... @@ -48,7 +54,7 @@ int main (int argc, char *argv[])
48 54 // Extract error message
49 55 std::string message = e.what();
50 56 std::string signal = "*** ERROR encounterted in the execution of"
51   - " ctlike. Run aborted ...";
  57 + " ctlike. Run aborted ...";
52 58  
53 59 // Write error in logger
54 60 application.log << signal << std::endl;
... ...