Konstancja Satalecka / ctools

Commit 13d498b918b48ba868c1e3d680f78bd4295138a2

Authored by knodlseder
1 parent 686648ba

generalize code and implement logger interface

Inline Side-by-side
Showing 2 changed files with 155 additions and 65 deletions
src/ctlike/ctlike.cpp
  View file @13d498b
... ... @@ -34,11 +34,14 @@
34 34 /***********************************************************************//**
35 35 * @brief Void constructor
36 36 ***************************************************************************/
37   -ctlike::ctlike(void) : GApplication(CTLIKE_NAME, CTLIKE_VERSION)
  37 +ctlike::ctlike(void) : GApplication(CTLIKE_NAME, CTLIKE_VERSION)
38 38 {
39 39 // Initialise members
40 40 init_members();
41 41  
  42 + // Write header into logger
  43 + log_header();
  44 +
42 45 // Return
43 46 return;
44 47 }
... ... @@ -56,6 +59,9 @@ ctlike::ctlike(int argc, char *argv[]) :
56 59 // Initialise members
57 60 init_members();
58 61  
  62 + // Write header into logger
  63 + log_header();
  64 +
59 65 // Return
60 66 return;
61 67 }
... ... @@ -88,75 +94,165 @@ ctlike::~ctlike(void)
88 94  
89 95 /***********************************************************************//**
90 96 * @brief Run maximum likelihood analysis
  97 + *
  98 + * The following analysis steps are performed:
  99 + * 1. Read the parameters (and write them into logger)
  100 + * 2. Load observation
  101 + * 3. Setup models for optimizing
  102 + * 4. Optimize model (and write result into logger)
  103 + * 5. Write results out as XML model
91 104 ***************************************************************************/
92   -int ctlike::run(void)
  105 +void ctlike::run(void)
93 106 {
94   - // Initialise return code
95   - int rc = 0;
  107 + // Get parameters
  108 + get_parameters();
  109 +
  110 + // Write parameters into logger
  111 + if (logTerse()) {
  112 + log_parameters();
  113 + log << std::endl;
  114 + }
  115 +
  116 + // Load observation
  117 + if (m_method == "UNBINNED")
  118 + unbinned(m_evfile);
  119 + else
  120 + binned(m_cntmap);
  121 +
  122 + // Setup models for optimizing.
  123 + m_models = GModels(m_srcmdl);
  124 + m_obs.models(m_models);
  125 +
  126 + // Optimize
  127 + optimize_lm();
  128 +
  129 + // Write results into logger
  130 + if (logTerse()) {
  131 + log << " Maximum log likelihood ....: " << m_logL << std::endl;
  132 + log << " Npred .....................: " << m_obs.npred() << std::endl;
  133 + log << m_models << std::endl << std::endl;
  134 + }
  135 +
  136 + // Write results out as XML model
  137 + if (toupper(m_outmdl) != "NONE")
  138 + m_models.save(m_outmdl);
  139 +
  140 + // Return
  141 + return;
  142 +}
96 143  
  144 +
  145 +/***********************************************************************//**
  146 + * @brief Get application parameters
  147 + *
  148 + * Get all task parameters from parameter file or (if required) by querying
  149 + * the user.
  150 + ***************************************************************************/
  151 +void ctlike::get_parameters(void)
  152 +{
97 153 // Get standard parameters
98 154 m_method = toupper(par("method")->value());
99 155 m_caldb = par("caldb")->value();
100 156 m_irf = par("irf")->value();
101 157 m_srcmdl = par("srcmdl")->value();
102 158 m_outmdl = par("outmdl")->value();
103   - m_models = GModels(m_srcmdl);
104 159  
105   - // Set fixed parameters
106   - m_opt = new GOptimizerLM;
  160 + // Get unbinned parameters ...
  161 + if (m_method == "UNBINNED") {
  162 + m_evfile = par("evfile")->value();
  163 + }
  164 +
  165 + // ... or get binned parameters
  166 + else if (m_method == "BINNED") {
  167 + m_cntmap = par("cntmap")->value();
  168 + }
  169 +
  170 + // Return
  171 + return;
  172 +}
  173 +
  174 +
  175 +/***********************************************************************//**
  176 + * @brief Optimize using Levenberg-Marquardt method
  177 + ***************************************************************************/
  178 +void ctlike::optimize_lm(void)
  179 +{
  180 + // Allocate optimizer. The logger is only passed to the optimizer
  181 + // constructor if optimizer logging is requested.
  182 + GOptimizerLM* opt = (logExplicit()) ? new GOptimizerLM(log)
  183 + : new GOptimizerLM();
  184 +
  185 + // Assign optimizer
  186 + m_opt = opt;
107 187  
108 188 // Set optimizer parameters
109   - ((GOptimizerLM*)m_opt)->max_iter(m_max_iter);
  189 + opt->max_iter(m_max_iter);
110 190  
111   - // Branch on method
112   - if (m_method == "UNBINNED")
113   - rc = unbinned();
114   - else
115   - rc = binned();
  191 + // Write Header for optimization and indent for optimizer logging
  192 + if (logExplicit()) {
  193 + log.header1("Maximum likelihood optimisation");
  194 + log.indent(1);
  195 + }
  196 +
  197 + // Perform LM optimization
  198 + m_obs.optimize(*m_opt);
116 199  
117   - // Write result as XML model
118   - if (toupper(m_outmdl) != "NONE")
119   - m_models.save(m_outmdl);
  200 + // Get models back
  201 + m_models = *(m_obs.models());
  202 +
  203 + // Store maximum log likelihood value
  204 + m_logL = -opt->value();
120 205  
121   - // Dump results
122   - std::cout << m_models << std::endl;
  206 + // Write optimization results
  207 + log.indent(0);
  208 + if (logTerse()) {
  209 + log << std::endl;
  210 + log.header1("Maximum likelihood optimisation results");
  211 + log << *opt << std::endl;
  212 + }
123 213  
124 214 // Free optimizer
125   - delete m_opt;
  215 + delete opt;
  216 +
  217 + // Reset optimizer pointer
  218 + m_opt = NULL;
126 219  
127 220 // Return
128   - return rc;
  221 + return;
129 222 }
130 223  
131 224  
132 225 /***********************************************************************//**
133   - * @brief Perform unbinned maximum likelihood analysis
  226 + * @brief Load unbinned observation
  227 + *
  228 + * @param[in] evfile Events file name.
  229 + *
  230 + * @todo ROI, energy range and time range information is hardwired. This is
  231 + * not needed in the future once this information is correctly extracted
  232 + * from the event file header.
134 233 ***************************************************************************/
135   -int ctlike::unbinned(void)
  234 +void ctlike::unbinned(const std::string& evfile)
136 235 {
137   - // Initialise return code
138   - int rc = 0;
139   -
140   - // Get parameters
141   - std::string evfile = par("evfile")->value();
142   -
143 236 // Declare observations
144 237 GCTAObservation run;
145 238  
146   - // Setup ROI covered by data
  239 + // DUMMY: Setup ROI covered by data. This is needed since actual test
  240 + // data do not contain any ROI informations.
147 241 GCTAInstDir instDir;
148 242 GCTARoi roi;
149 243 instDir.radec_deg(117.02, -33.35); // Adapt to file
150 244 roi.centre(instDir);
151 245 roi.radius(2.5);
152 246  
153   - // Setup energy range covered by data
  247 + // DUMMY: Setup energy range covered by data. This is needed since actual
  248 + // test data do not contain any energy range information
154 249 GEnergy emin;
155 250 GEnergy emax;
156 251 emin.TeV(0.02);
157 252 emax.TeV(100.0);
158 253  
159   - // Setup time range covered by data
  254 + // DUMMY: Setup time range covered by data. This is needed since actual
  255 + // test data do not contain any time range information
160 256 GTime tstart;
161 257 GTime tstop;
162 258 tstart.met(0.0);
... ... @@ -171,50 +267,32 @@ int ctlike::unbinned(void)
171 267 run.ebounds()->append(emin, emax);
172 268 m_obs.append(run);
173 269  
174   - // Setup models for optimizing
175   - m_obs.models(m_models);
176   -
177   - // Perform LM optimization
178   - m_obs.optimize(*m_opt);
179   -
180   - // Get models back
181   - m_models = *(m_obs.models());
182   -
183 270 // Return
184   - return rc;
  271 + return;
185 272 }
186 273  
187 274  
188 275 /***********************************************************************//**
189   - * @brief Perform binned maximum likelihood analysis
  276 + * @brief Load binned observation
  277 + *
  278 + * @param[in] cntmap Counts map file name.
190 279 ***************************************************************************/
191   -int ctlike::binned(void)
  280 +void ctlike::binned(const std::string& cntmap)
192 281 {
193   - // Initialise return code
194   - int rc = 0;
195   -
196   - // Get parameters
197   - std::string cntmap = par("cntmap")->value();
198   -
199 282 // Declare observations
200 283 GCTAObservation run;
201 284  
202 285 // Load binned CTA observation
203 286 run.load_binned(cntmap);
204   - run.response(m_irf, m_caldb);
205   - m_obs.append(run);
206 287  
207   - // Setup models for optimizing
208   - m_obs.models(m_models);
209   -
210   - // Perform LM optimization
211   - m_obs.optimize(*m_opt);
  288 + // Set response
  289 + run.response(m_irf, m_caldb);
212 290  
213   - // Get models back
214   - m_models = *(m_obs.models());
  291 + // Append observation to container
  292 + m_obs.append(run);
215 293  
216 294 // Return
217   - return rc;
  295 + return;
218 296 }
219 297  
220 298  
... ... @@ -236,10 +314,14 @@ void ctlike::init_members(void)
236 314 m_srcmdl.clear();
237 315 m_outmdl.clear();
238 316 m_models.clear();
239   - //m_obs.clear(); //!< NOT YET IMPLEMENTED
  317 + m_obs.clear();
240 318 m_max_iter = 1000;
  319 + m_logL = 0.0;
241 320 m_opt = NULL;
242 321  
  322 + // Set logger properties
  323 + log.date(true);
  324 +
243 325 // Return
244 326 return;
245 327 }
... ... @@ -251,7 +333,6 @@ void ctlike::init_members(void)
251 333 void ctlike::free_members(void)
252 334 {
253 335 // Free members
254   - delete m_opt;
255 336  
256 337 // Return
257 338 return;
... ...
src/ctlike/ctlike.hpp
  View file @13d498b
... ... @@ -42,9 +42,11 @@ public:
42 42 ~ctlike(void);
43 43  
44 44 // Methods
45   - int run(void);
46   - int unbinned(void);
47   - int binned(void);
  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);
48 50  
49 51 protected:
50 52 // Protected methods
... ... @@ -60,8 +62,15 @@ protected:
60 62 GModels m_models; //!< Source models
61 63 GObservations m_obs; //!< Observations
62 64  
63   - // Fixed parameters
  65 + // Unbinned likelihood parameters
  66 + std::string m_evfile; //!< Events file
  67 +
  68 + // Binned likelihood parameters
  69 + std::string m_cntmap; //!< Counts map
  70 +
  71 + // Internal parameters
64 72 int m_max_iter; //!< Maximum number of iterations
  73 + double m_logL; //!< Maximum log likelihood
65 74 GOptimizer* m_opt; //!< Optimizer
66 75 };
67 76  
... ...