Fabio Iocco / ctools

Commit 2aa78ba1527bfa6c047c4c422d744ca2923bac4e

Authored by Jürgen Knödlseder
1 parent a35f9317

Introduce ctlikelihood base class for likelihood ctools (#1828) 

All likelihood ctools now derive from the ctlikelihood base class. This avoids code duplication and provides all likelihood tools with a uniform interface.
Inline Side-by-side
Showing 26 changed files with 660 additions and 397 deletions
ChangeLog
  View file @2aa78ba
... ... @@ -3,6 +3,7 @@
3 3 * Version 1.2.0 released
4 4 ========================
5 5  
  6 + Introduce ctlikelihood base class (#1828)
6 7 Remove ctselect::time_reference() method
7 8 Add ctool and cscript pars() methods (#1838)
8 9 Fix ctbin::save() bug for empty observation container (#1835)
... ...
... ... @@ -24,7 +24,10 @@ all parameters in a dictionary (#1838).
24 24  
25 25 Support model
26 26 -------------
27   -None
  27 +The ctlikelihood base class was introduced from which all likelihood ctools
  28 +now derive (ctlike, cttsmap, cterror, ctulimit). This avoids code duplication
  29 +and provides all tools with a standard interface. All tools now have the
  30 +obs() and opt() method (#1828).
28 31  
29 32  
30 33 ctobssim - CTA observation simulation
... ... @@ -57,12 +60,16 @@ ctlike - Maximum likelihood fitting
57 60 The clear() method now calls GApplication::clear() to conserve the properties
58 61 of the tool (#1835).
59 62  
  63 +The class now derives from ctlikelihood (#1828).
  64 +
60 65  
61 66 cttsmap - Generation of Test Statistic map
62 67 ------------------------------------------
63 68 The clear() method now calls GApplication::clear() to conserve the properties
64 69 of the tool (#1835).
65 70  
  71 +The class now derives from ctlikelihood (#1828).
  72 +
66 73  
67 74 ctskymap - CTA sky mapping
68 75 --------------------------
... ... @@ -134,12 +141,16 @@ of the tool (#1835).
134 141 The tool now also supports spectral models with the new naming conventions, and
135 142 in particular the "PhotonFlux" and "EnergyFlux" names (#1832).
136 143  
  144 +The class now derives from ctlikelihood (#1828).
  145 +
137 146  
138 147 cterror - Likelihood profile error computation
139 148 ----------------------------------------------
140 149 The clear() method now calls GApplication::clear() to conserve the properties
141 150 of the tool (#1835).
142 151  
  152 +The class now derives from ctlikelihood (#1828).
  153 +
143 154  
144 155 cscripts
145 156 --------
... ...
pyext/Makefile.am
  View file @2aa78ba
... ... @@ -168,8 +168,8 @@ CLEANFILES = $(CLEAN_FILES) $(DEPFILES)
168 168 # Files to be removed in build directory
169 169 clean-local:
170 170 rm -rf build
171   - rm -rf $(builddir)/ctools/*.py
172   - rm -rf $(builddir)/ctools/*.pyc
  171 + rm -rf $(builddir)/ctools/c*.py
  172 + rm -rf $(builddir)/ctools/c*.pyc
173 173 rm -rf $(builddir)/ctools/*_wrap.cpp
174 174 rm -rf $(builddir)/ctools/_tools.so
175 175 rm -rf $(builddir)/ctools/__pycache__
... ...
pyext/cterror.i
  View file @2aa78ba
1 1 /***************************************************************************
2 2 * cterror - Parameter error calculation tool *
3 3 * ----------------------------------------------------------------------- *
4   - * copyright (C) 2015 by Florent Forest *
  4 + * copyright (C) 2015-2016 by Florent Forest *
5 5 * ----------------------------------------------------------------------- *
6 6 * *
7 7 * This program is free software: you can redistribute it and/or modify *
... ... @@ -34,9 +34,8 @@
34 34 * @class cterror
35 35 *
36 36 * @brief Parameter error calculation tool
37   - *
38 37 ***************************************************************************/
39   -class cterror : public ctool {
  38 +class cterror : public ctlikelihood {
40 39  
41 40 public:
42 41 // Constructors and destructors
... ... @@ -47,11 +46,9 @@ public:
47 46 virtual ~cterror(void);
48 47  
49 48 // Methods
50   - void clear(void);
51   - void run(void);
52   - void save(void);
53   - const GObservations& obs(void) const;
54   - const GOptimizer* opt(void) const;
  49 + void clear(void);
  50 + void run(void);
  51 + void save(void);
55 52 };
56 53  
57 54  
... ...
pyext/ctlike.i
  View file @2aa78ba
1 1 /***************************************************************************
2 2 * ctlike - Maximum likelihood fitting tool *
3 3 * ----------------------------------------------------------------------- *
4   - * copyright (C) 2010-2014 by Jurgen Knodlseder *
  4 + * copyright (C) 2010-2016 by Jurgen Knodlseder *
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 * @brief Maximum likelihood fitting tool
36 36 ***************************************************************************/
37   -class ctlike : public ctool {
  37 +class ctlike : public ctlikelihood {
38 38 public:
39 39 // Constructors and destructors
40 40 ctlike(void);
... ... @@ -44,11 +44,9 @@ public:
44 44 virtual ~ctlike(void);
45 45  
46 46 // Methods
47   - void clear(void);
48   - void run(void);
49   - void save(void);
50   - const GObservations& obs(void) const;
51   - const GOptimizer* opt(void) const;
  47 + void clear(void);
  48 + void run(void);
  49 + void save(void);
52 50 };
53 51  
54 52  
... ...
pyext/ctlikelihood.i 0 → 100644
  View file @2aa78ba
  1 +/***************************************************************************
  2 + * ctlikelihood - Base class for likelihood tools *
  3 + * ----------------------------------------------------------------------- *
  4 + * copyright (C) 2016 by Juergen Knoedlseder *
  5 + * ----------------------------------------------------------------------- *
  6 + * *
  7 + * This program is free software: you can redistribute it and/or modify *
  8 + * it under the terms of the GNU General Public License as published by *
  9 + * the Free Software Foundation, either version 3 of the License, or *
  10 + * (at your option) any later version. *
  11 + * *
  12 + * This program is distributed in the hope that it will be useful, *
  13 + * but WITHOUT ANY WARRANTY; without even the implied warranty of *
  14 + * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the *
  15 + * GNU General Public License for more details. *
  16 + * *
  17 + * You should have received a copy of the GNU General Public License *
  18 + * along with this program. If not, see <http://www.gnu.org/licenses/>. *
  19 + * *
  20 + ***************************************************************************/
  21 +/**
  22 + * @file ctlikelihood.i
  23 + * @brief Likelihood tool base class interface definition
  24 + * @author Juergen Knoedlseder
  25 + */
  26 +
  27 +%{
  28 +/* Put headers and other declarations here that are needed for compilation */
  29 +#include "ctlikelihood.hpp"
  30 +%}
  31 +
  32 +
  33 +/***********************************************************************//**
  34 + * @class ctlikelihood
  35 + *
  36 + * @brief Base class for likelihood tools
  37 + ***************************************************************************/
  38 +class ctlikelihood : public ctool {
  39 +
  40 +public:
  41 + // Constructors and destructors
  42 + ctlikelihood(void);
  43 + ctlikelihood(const std::string& name, const std::string& version);
  44 + ctlikelihood(const std::string& name, const std::string& version,
  45 + const GObservations& obs);
  46 + ctlikelihood(const std::string& name, const std::string& version,
  47 + int argc, char* argv[]);
  48 + ctlikelihood(const ctlikelihood& app);
  49 + virtual ~ctlikelihood(void);
  50 +
  51 + // Pure virtual methods
  52 + virtual void clear(void) = 0;
  53 + virtual void run(void) = 0;
  54 + virtual void save(void) = 0;
  55 +
  56 + // Methods
  57 + const GObservations& obs(void) const;
  58 + const GOptimizer* opt(void) const;
  59 +};
  60 +
  61 +
  62 +/***********************************************************************//**
  63 + * @brief Likelihood tool Python extensions
  64 + ***************************************************************************/
  65 +%extend ctlikelihood {
  66 +}
... ...
pyext/ctool.i
  View file @2aa78ba
... ... @@ -62,7 +62,7 @@ public:
62 62 /***********************************************************************//**
63 63 * @class ctool
64 64 *
65   - * @brief ctool base class
  65 + * @brief Base class for ctools
66 66 ***************************************************************************/
67 67 class ctool : public GApplication {
68 68 public:
... ...
pyext/ctools/tools.i
  View file @2aa78ba
... ... @@ -44,8 +44,11 @@
44 44 %import(module="gammalib.base") "GBase.i";
45 45 %import(module="gammalib.app") "GApplication.i";
46 46  
47   -/* __ ctools _____________________________________________________________ */
  47 +/* __ ctools base classes ________________________________________________ */
48 48 %include "ctool.i"
  49 +%include "ctlikelihood.i"
  50 +
  51 +/* __ ctools _____________________________________________________________ */
49 52 %include "ctbin.i"
50 53 %include "ctobssim.i"
51 54 %include "ctlike.i"
... ...
pyext/cttsmap.i
  View file @2aa78ba
... ... @@ -34,7 +34,8 @@
34 34 *
35 35 * @brief TS map calculation tool
36 36 ***************************************************************************/
37   -class cttsmap : public ctool {
  37 +class cttsmap : public ctlikelihood {
  38 +
38 39 public:
39 40 // Constructors and destructors
40 41 cttsmap(void);
... ... @@ -44,12 +45,11 @@ public:
44 45 virtual ~cttsmap(void);
45 46  
46 47 // Methods
47   - void clear(void);
48   - void run(void);
49   - void save(void);
50   - void publish(const std::string& name = "");
51   - const GObservations& obs(void) const;
52   - const GSkyMap& tsmap(void) const;
  48 + void clear(void);
  49 + void run(void);
  50 + void save(void);
  51 + void publish(const std::string& name = "");
  52 + const GSkyMap& tsmap(void) const;
53 53 };
54 54  
55 55  
... ...
pyext/ctulimit.i
  View file @2aa78ba
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 * @brief Upper limit calculation tool
36 36 ***************************************************************************/
37   -class ctulimit : public ctool {
  37 +class ctulimit : public ctlikelihood {
38 38 public:
39 39 // Constructors and destructors
40 40 ctulimit(void);
... ... @@ -44,14 +44,12 @@ public:
44 44 virtual ~ctulimit(void);
45 45  
46 46 // Methods
47   - void clear(void);
48   - void run(void);
49   - void save(void);
50   - const GObservations& obs(void) const;
51   - const GOptimizer* opt(void) const;
52   - const double& diff_ulimit(void) const;
53   - const double& flux_ulimit(void) const;
54   - const double& eflux_ulimit(void) const;
  47 + void clear(void);
  48 + void run(void);
  49 + void save(void);
  50 + const double& diff_ulimit(void) const;
  51 + const double& flux_ulimit(void) const;
  52 + const double& eflux_ulimit(void) const;
55 53 };
56 54  
57 55  
... ...
src/cterror/cterror.cpp
  View file @2aa78ba
... ... @@ -35,7 +35,6 @@
35 35  
36 36 /* __ Method name definitions ____________________________________________ */
37 37 #define G_ERR_BISECTION "cterror::error_bisection(double&, double&)"
38   -#define G_EVALUATE "cterror::evaluate(double&)"
39 38  
40 39 /* __ Debug definitions __________________________________________________ */
41 40  
... ... @@ -53,7 +52,7 @@
53 52 *
54 53 * Constructs an empty cterror tool.
55 54 ***************************************************************************/
56   -cterror::cterror(void) : ctool(CTERROR_NAME, CTERROR_VERSION)
  55 +cterror::cterror(void) : ctlikelihood(CTERROR_NAME, CTERROR_VERSION)
57 56 {
58 57 // Initialise members
59 58 init_members();
... ... @@ -68,24 +67,19 @@ cterror::cterror(void) : ctool(CTERROR_NAME, CTERROR_VERSION)
68 67 *
69 68 * param[in] obs Observation container.
70 69 *
71   - * Constructs an instance of the cterror tool that will operate on the
72   - * provided observation container.
  70 + * Constructs cterror tool from an observations container.
73 71 ***************************************************************************/
74 72 cterror::cterror(const GObservations& obs) :
75   - ctool(CTERROR_NAME, CTERROR_VERSION)
  73 + ctlikelihood(CTERROR_NAME, CTERROR_VERSION, obs)
76 74 {
77 75 // Initialise members
78 76 init_members();
79 77  
80   - // Set observations
81   - m_obs = obs;
82   -
83 78 // Return
84 79 return;
85 80 }
86 81  
87 82  
88   -
89 83 /***********************************************************************//**
90 84 * @brief Command line constructor
91 85 *
... ... @@ -96,7 +90,7 @@ cterror::cterror(const GObservations&amp; obs) :
96 90 * parameters that are provided as command line arguments.
97 91 ***************************************************************************/
98 92 cterror::cterror(int argc, char *argv[]) :
99   - ctool(CTERROR_NAME, CTERROR_VERSION, argc, argv)
  93 + ctlikelihood(CTERROR_NAME, CTERROR_VERSION, argc, argv)
100 94 {
101 95 // Initialise members
102 96 init_members();
... ... @@ -114,7 +108,7 @@ cterror::cterror(int argc, char *argv[]) :
114 108 * Constructs an instance of the cterror tool by copying information from
115 109 * another ctulimit tool.
116 110 ***************************************************************************/
117   -cterror::cterror(const cterror& app) : ctool(app)
  111 +cterror::cterror(const cterror& app) : ctlikelihood(app)
118 112 {
119 113 // Initialise members
120 114 init_members();
... ... @@ -162,7 +156,7 @@ cterror&amp; cterror::operator=(const cterror&amp; app)
162 156 if (this != &app) {
163 157  
164 158 // Copy base class members
165   - this->ctool::operator=(app);
  159 + this->ctlikelihood::operator=(app);
166 160  
167 161 // Free members
168 162 free_members();
... ... @@ -195,6 +189,7 @@ void cterror::clear(void)
195 189 {
196 190 // Free members
197 191 free_members();
  192 + this->ctlikelihood::free_members();
198 193 this->ctool::free_members();
199 194  
200 195 // Clear base class (needed to conserve tool name and version)
... ... @@ -202,6 +197,7 @@ void cterror::clear(void)
202 197  
203 198 // Initialise members
204 199 this->ctool::init_members();
  200 + this->ctlikelihood::init_members();
205 201 init_members();
206 202  
207 203 // Write header into logger
... ... @@ -384,19 +380,17 @@ void cterror::init_members(void)
384 380 // Initialise user parameters
385 381 m_srcname.clear();
386 382 m_outmodel.clear();
387   - m_confidence = 0.68;
388   - m_tol = 1.0e-3;
389   - m_max_iter = 50;
390   - m_value = 0.0;
391   - m_chatter = static_cast<GChatter>(2);
  383 + m_confidence = 0.68;
  384 + m_tol = 1.0e-3;
  385 + m_max_iter = 50;
  386 + m_apply_edisp = false;
  387 + m_chatter = static_cast<GChatter>(2);
392 388  
393 389 // Initialise protected members
394   - m_obs.clear();
395   - m_opt.clear();
  390 + m_value = 0.0;
396 391 m_dlogL = 0.0;
397 392 m_best_logL = 0.0;
398 393 m_model_par = NULL;
399   - m_apply_edisp = false;
400 394  
401 395 // Set optimizer parameters
402 396 m_opt.max_iter(m_max_iter);
... ... @@ -424,10 +418,9 @@ void cterror::copy_members(const cterror&amp; app)
424 418 m_chatter = app.m_chatter;
425 419  
426 420 // Copy protected members
427   - m_obs = app.m_obs;
  421 + m_value = app.m_value;
428 422 m_dlogL = app.m_dlogL;
429 423 m_best_logL = app.m_best_logL;
430   - m_opt = app.m_opt;
431 424 m_model_par = NULL;
432 425  
433 426 // Return
... ... @@ -548,7 +541,7 @@ double cterror::error_bisection(const double&amp; min, const double&amp; max)
548 541 mid = (wrk_min + wrk_max) / 2.0;
549 542  
550 543 // Calculate function value
551   - double eval_mid = evaluate(mid);
  544 + double eval_mid = evaluate(*m_model_par, mid) - (m_best_logL + m_dlogL);
552 545  
553 546 // Write interval into logger
554 547 log_value(EXPLICIT, " Iteration "+gammalib::str(iter),
... ... @@ -597,53 +590,3 @@ double cterror::error_bisection(const double&amp; min, const double&amp; max)
597 590 return mid;
598 591  
599 592 }
600   -
601   -
602   -/***********************************************************************//**
603   - * @brief Evaluates the log-likelihood
604   - *
605   - * @param[in] value Parameter factor value
606   - * @return Log-likelihood value
607   - *
608   - * Evaluates the log-likelihood as a function of the parameter of interest.
609   - ***************************************************************************/
610   -double cterror::evaluate(const double& value)
611   -{
612   - // Initialise log-likelihood value
613   - double logL = 0.0;
614   -
615   - // Check if given parameter is within boundaries
616   - if (value > m_model_par->factor_min() && value < m_model_par->factor_max()) {
617   -
618   - // Change parameter factor
619   - m_model_par->factor_value(value);
620   -
621   - // Fix parameter
622   - m_model_par->fix();
623   -
624   - // Re-optimize
625   - m_obs.optimize(m_opt);
626   -
627   - // Retrieve likelihood
628   - logL = m_obs.logL();
629   -
630   - // Free parameter
631   - m_model_par->free();
632   -
633   - } // endif: value was inside allowed range
634   -
635   - // ... otherwise signal that the parameter went outside the boundaries
636   - else {
637   - std::string msg = "Value of parameter \""+m_model_par->name()+"\" "
638   - "outside of validity range requested. To omit "
639   - "this error please enlarge the parameter range "
640   - "in the model XML file.";
641   - throw GException::invalid_value(G_EVALUATE, msg);
642   - }
643   -
644   - // Compute function value
645   - double logL_difference = logL - m_best_logL - m_dlogL;
646   -
647   - // Return
648   - return logL_difference;
649   -}
... ...
src/cterror/cterror.hpp
  View file @2aa78ba
... ... @@ -28,10 +28,7 @@
28 28 #define CTERROR_HPP
29 29  
30 30 /* __ Includes ___________________________________________________________ */
31   -#include <string>
32   -#include "GammaLib.hpp"
33   -#include "GCTALib.hpp"
34   -#include "ctool.hpp"
  31 +#include "ctlikelihood.hpp"
35 32  
36 33 /* __Definitions _________________________________________________________ */
37 34 #define CTERROR_NAME "cterror"
... ... @@ -43,7 +40,7 @@
43 40 *
44 41 * @brief Parameter error calculation tool
45 42 ***************************************************************************/
46   -class cterror : public ctool {
  43 +class cterror : public ctlikelihood {
47 44  
48 45 public:
49 46 // Constructors and destructors
... ... @@ -57,11 +54,9 @@ public:
57 54 cterror& operator=(const cterror& app);
58 55  
59 56 // Methods
60   - void clear(void);
61   - void run(void);
62   - void save(void);
63   - const GObservations& obs(void) const;
64   - const GOptimizer* opt(void) const;
  57 + void clear(void);
  58 + void run(void);
  59 + void save(void);
65 60  
66 61 protected:
67 62 // Protected methods
... ... @@ -69,7 +64,6 @@ protected:
69 64 void copy_members(const cterror& app);
70 65 void free_members(void);
71 66 void get_parameters(void);
72   - double evaluate(const double& value);
73 67 double error_bisection(const double& min, const double& max);
74 68  
75 69 // User parameters
... ... @@ -78,40 +72,14 @@ protected:
78 72 double m_confidence; //!< Confidence level
79 73 double m_tol; //!< Tolerance for limit determination
80 74 int m_max_iter; //!< Maximum number of iterations
81   - double m_value; //!< Parameter value
82 75 bool m_apply_edisp; //!< Apply energy dispersion?
83 76 GChatter m_chatter; //!< Chattiness
84 77  
85 78 // Protected members
86   - GObservations m_obs; //!< Observation container
  79 + double m_value; //!< Parameter value
87 80 double m_dlogL; //!< Likelihood difference for upper limit computation
88 81 GModelPar* m_model_par; //!< Pointer to model parameter
89 82 double m_best_logL; //!< Best fit log likelihood of given model
90   - GOptimizerLM m_opt; //!< Optimizer
91 83 };
92 84  
93   -
94   -/***********************************************************************//**
95   - * @brief Return observation container
96   - *
97   - * @return Reference to observation container.
98   - ***************************************************************************/
99   -inline
100   -const GObservations& cterror::obs(void) const
101   -{
102   - return m_obs;
103   -}
104   -
105   -
106   -/***********************************************************************//**
107   - * @brief Return optimizer
108   - *
109   - * @return Pointer to optimizer.
110   - ***************************************************************************/
111   -inline
112   -const GOptimizer* cterror::opt(void) const
113   -{
114   - return &m_opt;
115   -}
116   -
117 85 #endif /* CTERROR_HPP */
... ...
src/ctlike/ctlike.cpp
  View file @2aa78ba
... ... @@ -48,7 +48,7 @@
48 48 /***********************************************************************//**
49 49 * @brief Void constructor
50 50 ***************************************************************************/
51   -ctlike::ctlike(void) : ctool(CTLIKE_NAME, CTLIKE_VERSION)
  51 +ctlike::ctlike(void) : ctlikelihood(CTLIKE_NAME, CTLIKE_VERSION)
52 52 {
53 53 // Initialise members
54 54 init_members();
... ... @@ -63,17 +63,14 @@ ctlike::ctlike(void) : ctool(CTLIKE_NAME, CTLIKE_VERSION)
63 63 *
64 64 * param[in] obs Observation container.
65 65 *
66   - * This method creates an instance of the class by copying an existing
67   - * observations container.
  66 + * Constructs ctlike tool from an observations container.
68 67 ***************************************************************************/
69   -ctlike::ctlike(const GObservations& obs) : ctool(CTLIKE_NAME, CTLIKE_VERSION)
  68 +ctlike::ctlike(const GObservations& obs) :
  69 + ctlikelihood(CTLIKE_NAME, CTLIKE_VERSION, obs)
70 70 {
71 71 // Initialise members
72 72 init_members();
73 73  
74   - // Set observations
75   - m_obs = obs;
76   -
77 74 // Return
78 75 return;
79 76 }
... ... @@ -84,9 +81,12 @@ ctlike::ctlike(const GObservations&amp; obs) : ctool(CTLIKE_NAME, CTLIKE_VERSION)
84 81 *
85 82 * @param[in] argc Number of arguments in command line.
86 83 * @param[in] argv Array of command line arguments.
  84 + *
  85 + * Constructs an instance of the ctlike tool that will parse user parameters
  86 + * that are provided as command line arguments.
87 87 ***************************************************************************/
88 88 ctlike::ctlike(int argc, char *argv[]) :
89   - ctool(CTLIKE_NAME, CTLIKE_VERSION, argc, argv)
  89 + ctlikelihood(CTLIKE_NAME, CTLIKE_VERSION, argc, argv)
90 90 {
91 91 // Initialise members
92 92 init_members();
... ... @@ -101,7 +101,7 @@ ctlike::ctlike(int argc, char *argv[]) :
101 101 *
102 102 * @param[in] app Application.
103 103 ***************************************************************************/
104   -ctlike::ctlike(const ctlike& app) : ctool(app)
  104 +ctlike::ctlike(const ctlike& app) : ctlikelihood(app)
105 105 {
106 106 // Initialise members
107 107 init_members();
... ... @@ -145,7 +145,7 @@ ctlike&amp; ctlike::operator=(const ctlike&amp; app)
145 145 if (this != &app) {
146 146  
147 147 // Copy base class members
148   - this->ctool::operator=(app);
  148 + this->ctlikelihood::operator=(app);
149 149  
150 150 // Free members
151 151 free_members();
... ... @@ -178,6 +178,7 @@ void ctlike::clear(void)
178 178 {
179 179 // Free members
180 180 free_members();
  181 + this->ctlikelihood::free_members();
181 182 this->ctool::free_members();
182 183  
183 184 // Clear base class (needed to conserve tool name and version)
... ... @@ -185,6 +186,7 @@ void ctlike::clear(void)
185 186  
186 187 // Initialise members
187 188 this->ctool::init_members();
  189 + this->ctlikelihood::init_members();
188 190 init_members();
189 191  
190 192 // Write header into logger
... ... @@ -291,7 +293,7 @@ void ctlike::run(void)
291 293  
292 294 // Write results into logger
293 295 log_header1(NORMAL, "Maximum likelihood optimisation results");
294   - log_string(NORMAL, m_opt->print(m_chatter));
  296 + log_string(NORMAL, m_opt.print(m_chatter));
295 297 log_value(NORMAL, "Maximum log likelihood", gammalib::str(m_logL,3));
296 298 log_value(NORMAL, "Observed events (Nobs)", gammalib::str(num_events,3));
297 299 log_value(NORMAL, "Predicted events (Npred)", gammalib::str(npred,3)+
... ... @@ -398,10 +400,10 @@ void ctlike::get_parameters(void)
398 400  
399 401 // Set optimizer logger
400 402 if (logNormal()) {
401   - static_cast<GOptimizerLM*>(m_opt)->logger(&log);
  403 + static_cast<GOptimizerLM*>(&m_opt)->logger(&log);
402 404 }
403 405 else {
404   - static_cast<GOptimizerLM*>(m_opt)->logger(NULL);
  406 + static_cast<GOptimizerLM*>(&m_opt)->logger(NULL);
405 407 }
406 408  
407 409 // Write parameters into logger
... ... @@ -442,7 +444,7 @@ void ctlike::optimize_lm(void)
442 444 }
443 445  
444 446 // Perform LM optimization
445   - m_obs.optimize(*m_opt);
  447 + m_obs.optimize(m_opt);
446 448  
447 449 // Optionally refit
448 450 if (m_refit) {
... ... @@ -453,7 +455,7 @@ void ctlike::optimize_lm(void)
453 455 log.indent(1);
454 456  
455 457 // Optimise again
456   - m_obs.optimize(*m_opt);
  458 + m_obs.optimize(m_opt);
457 459  
458 460 }
459 461  
... ... @@ -464,10 +466,10 @@ void ctlike::optimize_lm(void)
464 466 (m_obs.function()).curvature())->print());
465 467  
466 468 // Compute errors
467   - m_obs.errors(*m_opt);
  469 + m_obs.errors(m_opt);
468 470  
469 471 // Store maximum log likelihood value
470   - m_logL = -(m_opt->value());
  472 + m_logL = -(m_opt.value());
471 473  
472 474 // Remove indent
473 475 log.indent(0);
... ... @@ -490,7 +492,7 @@ double ctlike::reoptimize_lm(void)
490 492 log.indent(1);
491 493  
492 494 // Create a clone of the optimizer for the re-optimisation
493   - GOptimizer* opt = m_opt->clone();
  495 + GOptimizer* opt = m_opt.clone();
494 496  
495 497 // Perform LM optimization
496 498 m_obs.optimize(*opt);
... ... @@ -526,12 +528,10 @@ void ctlike::init_members(void)
526 528 {
527 529 // Initialise members
528 530 m_outmodel.clear();
529   - m_obs.clear();
530 531 m_refit = false;
531 532 m_max_iter = 100; // Set maximum number of iterations
532 533 m_max_stall = 10; // Set maximum number of stalls
533 534 m_logL = 0.0;
534   - m_opt = NULL;
535 535 m_apply_edisp = false;
536 536 m_fix_spat_for_ts = false;
537 537 m_chatter = static_cast<GChatter>(2);
... ... @@ -539,15 +539,9 @@ void ctlike::init_members(void)
539 539 // Set logger properties
540 540 log.date(true);
541 541  
542   - // Allocate LM optimizer
543   - GOptimizerLM* opt = new GOptimizerLM();
544   -
545 542 // Set optimizer parameters
546   - opt->max_iter(m_max_iter);
547   - opt->max_stalls(m_max_stall);
548   -
549   - // Set optimizer pointer
550   - m_opt = opt;
  543 + m_opt.max_iter(m_max_iter);
  544 + m_opt.max_stalls(m_max_stall);
551 545  
552 546 // Return
553 547 return;
... ... @@ -564,11 +558,9 @@ void ctlike::copy_members(const ctlike&amp; app)
564 558 // Copy attributes
565 559 m_refit = app.m_refit;
566 560 m_outmodel = app.m_outmodel;
567   - m_obs = app.m_obs;
568 561 m_max_iter = app.m_max_iter;
569 562 m_max_stall = app.m_max_stall;
570 563 m_logL = app.m_logL;
571   - m_opt = app.m_opt->clone();
572 564 m_apply_edisp = app.m_apply_edisp;
573 565 m_fix_spat_for_ts = app.m_fix_spat_for_ts;
574 566 m_chatter = app.m_chatter;
... ... @@ -583,12 +575,6 @@ void ctlike::copy_members(const ctlike&amp; app)
583 575 ***************************************************************************/
584 576 void ctlike::free_members(void)
585 577 {
586   - // Free members
587   - if (m_opt != NULL) delete m_opt;
588   -
589   - // Mark pointers as free
590   - m_opt = NULL;
591   -
592 578 // Return
593 579 return;
594 580 }
... ...
src/ctlike/ctlike.hpp
  View file @2aa78ba
... ... @@ -28,14 +28,11 @@
28 28 #define CTLIKE_HPP
29 29  
30 30 /* __ Includes ___________________________________________________________ */
31   -#include <string>
32   -#include "GammaLib.hpp"
33   -#include "GCTALib.hpp"
34   -#include "ctool.hpp"
  31 +#include "ctlikelihood.hpp"
35 32  
36 33 /* __Definitions _________________________________________________________ */
37 34 #define CTLIKE_NAME "ctlike"
38   -#define CTLIKE_VERSION "1.1.0"
  35 +#define CTLIKE_VERSION "1.2.0"
39 36  
40 37  
41 38 /***********************************************************************//**
... ... @@ -43,7 +40,7 @@
43 40 *
44 41 * @brief Maximum likelihood fitting tool
45 42 ***************************************************************************/
46   -class ctlike : public ctool {
  43 +class ctlike : public ctlikelihood {
47 44  
48 45 public:
49 46 // Constructors and destructors
... ... @@ -57,11 +54,9 @@ public:
57 54 ctlike& operator=(const ctlike& app);
58 55  
59 56 // Methods
60   - void clear(void);
61   - void run(void);
62   - void save(void);
63   - const GObservations& obs(void) const;
64   - const GOptimizer* opt(void) const;
  57 + void clear(void);
  58 + void run(void);
  59 + void save(void);
65 60  
66 61 protected:
67 62 // Protected methods
... ... @@ -80,35 +75,9 @@ protected:
80 75 GChatter m_chatter; //!< Chattiness
81 76  
82 77 // Members
83   - GObservations m_obs; //!< Observations
84 78 int m_max_iter; //!< Maximum number of iterations
85 79 int m_max_stall; //!< Maximum number of stalls
86 80 double m_logL; //!< Maximum log likelihood
87   - GOptimizer* m_opt; //!< Optimizer
88 81 };
89 82  
90   -
91   -/***********************************************************************//**
92   - * @brief Return observation container
93   - *
94   - * @return Reference to observation container
95   - ***************************************************************************/
96   -inline
97   -const GObservations& ctlike::obs(void) const
98   -{
99   - return m_obs;
100   -}
101   -
102   -
103   -/***********************************************************************//**
104   - * @brief Return optimizer
105   - *
106   - * @return Pointer to optimizer
107   - ***************************************************************************/
108   -inline
109   -const GOptimizer* ctlike::opt(void) const
110   -{
111   - return m_opt;
112   -}
113   -
114 83 #endif /* CTLIKE_HPP */
... ...
src/ctmodel/ctmodel.cpp
  View file @2aa78ba
... ... @@ -598,7 +598,7 @@ void ctmodel::get_parameters(void)
598 598 * "expcube", "psfcube", and "bkgcube" parameters). In that case, it appended
599 599 * a dummy event cube to the observation.
600 600 *
601   - * @TODO Support stacked energy dispersion
  601 + * @todo Support stacked energy dispersion
602 602 ***************************************************************************/
603 603 void ctmodel::get_obs(void)
604 604 {
... ...
src/ctobssim/ctobssim.cpp
  View file @2aa78ba
... ... @@ -1039,10 +1039,11 @@ void ctobssim::simulate_interval(GCTAObservation* obs,
1039 1039 * @brief Simulate source events for a time slice
1040 1040 *
1041 1041 * @param[in] obs Pointer on CTA observation.
  1042 + * @param[in] rsp Pointer on CTA response.
1042 1043 * @param[in,out] events Pointer on CTA event list.
1043   - * @param[in] models Model list.
1044   - * @param[in] tmin Start time.
1045   - * @param[in] tmax Stop time.
  1044 + * @param[in] model Model.
  1045 + * @param[in] tstart Start time.
  1046 + * @param[in] tstop Stop time.
1046 1047 * @param[in] etrue_min Minimum true energy.
1047 1048 * @param[in] etrue_max Maximum true energy.
1048 1049 * @param[in] ereco_min Minimum reconstructed energy.
... ...
src/cttsmap/cttsmap.cpp
  View file @2aa78ba
... ... @@ -48,8 +48,10 @@
48 48  
49 49 /***********************************************************************//**
50 50 * @brief Void constructor
  51 + *
  52 + * Constructs an empty cttsmap tool.
51 53 ***************************************************************************/
52   -cttsmap::cttsmap(void) : ctool(CTTSMAP_NAME, CTTSMAP_VERSION)
  54 +cttsmap::cttsmap(void) : ctlikelihood(CTTSMAP_NAME, CTTSMAP_VERSION)
53 55 {
54 56 // Initialise members
55 57 init_members();
... ... @@ -64,32 +66,30 @@ cttsmap::cttsmap(void) : ctool(CTTSMAP_NAME, CTTSMAP_VERSION)
64 66 *
65 67 * param[in] obs Observation container.
66 68 *
67   - * This method creates an instance of the class by copying an existing
68   - * observations container.
  69 + * Constructs cttsmap tool from an observations container.
69 70 ***************************************************************************/
70 71 cttsmap::cttsmap(const GObservations& obs) :
71   - ctool(CTTSMAP_NAME, CTTSMAP_VERSION)
  72 + ctlikelihood(CTTSMAP_NAME, CTTSMAP_VERSION, obs)
72 73 {
73 74 // Initialise members
74 75 init_members();
75 76  
76   - // Set observations
77   - m_obs = obs;
78   -
79 77 // Return
80 78 return;
81 79 }
82 80  
83 81  
84   -
85 82 /***********************************************************************//**
86 83 * @brief Command line constructor
87 84 *
88 85 * @param[in] argc Number of arguments in command line.
89 86 * @param[in] argv Array of command line arguments.
  87 + *
  88 + * Constructs an instance of the cttsmap tool that will parse user parameters
  89 + * that are provided as command line arguments.
90 90 ***************************************************************************/
91 91 cttsmap::cttsmap(int argc, char *argv[]) :
92   - ctool(CTTSMAP_NAME, CTTSMAP_VERSION, argc, argv)
  92 + ctlikelihood(CTTSMAP_NAME, CTTSMAP_VERSION, argc, argv)
93 93 {
94 94 // Initialise members
95 95 init_members();
... ... @@ -104,7 +104,7 @@ cttsmap::cttsmap(int argc, char *argv[]) :
104 104 *
105 105 * @param[in] app Application.
106 106 ***************************************************************************/
107   -cttsmap::cttsmap(const cttsmap& app) : ctool(app)
  107 +cttsmap::cttsmap(const cttsmap& app) : ctlikelihood(app)
108 108 {
109 109 // Initialise members
110 110 init_members();
... ... @@ -148,7 +148,7 @@ cttsmap&amp; cttsmap::operator=(const cttsmap&amp; app)
148 148 if (this != &app) {
149 149  
150 150 // Copy base class members
151   - this->ctool::operator=(app);
  151 + this->ctlikelihood::operator=(app);
152 152  
153 153 // Free members
154 154 free_members();
... ... @@ -181,6 +181,7 @@ void cttsmap::clear(void)
181 181 {
182 182 // Free members
183 183 free_members();
  184 + this->ctlikelihood::free_members();
184 185 this->ctool::free_members();
185 186  
186 187 // Clear base class (needed to conserve tool name and version)
... ... @@ -188,6 +189,7 @@ void cttsmap::clear(void)
188 189  
189 190 // Initialise members
190 191 this->ctool::init_members();
  192 + this->ctlikelihood::init_members();
191 193 init_members();
192 194  
193 195 // Write header into logger
... ... @@ -232,9 +234,13 @@ void cttsmap::run(void)
232 234 int binmin = (m_binmin == -1) ? 0 : m_binmin;
233 235 int binmax = (m_binmax == -1) ? m_tsmap.npix() : m_binmax;
234 236  
235   - // Initialise optimizer
236   - GOptimizerLM* opt = (logExplicit()) ? new GOptimizerLM(&log)
237   - : new GOptimizerLM();
  237 + // Set optimizer logger
  238 + if (logExplicit()) {
  239 + static_cast<GOptimizerLM*>(&m_opt)->logger(&log);
  240 + }
  241 + else {
  242 + static_cast<GOptimizerLM*>(&m_opt)->logger(NULL);
  243 + }
238 244  
239 245 // Store initial models
240 246 GModels models_orig = m_obs.models();
... ... @@ -253,8 +259,8 @@ void cttsmap::run(void)
253 259  
254 260 // Compute likelihood without the test source
255 261 m_obs.models(models);
256   - m_obs.optimize(*opt);
257   - m_logL0 = -(opt->value());
  262 + m_obs.optimize(m_opt);
  263 + m_logL0 = -(m_opt.value());
258 264  
259 265 }
260 266  
... ... @@ -280,10 +286,10 @@ void cttsmap::run(void)
280 286 m_obs.models(models);
281 287  
282 288 // Optimize observation container
283   - m_obs.optimize(*opt);
  289 + m_obs.optimize(m_opt);
284 290  
285 291 // Retrieve the Likelihood value
286   - double logL1 = -(opt->value());
  292 + double logL1 = -(m_opt.value());
287 293  
288 294 // Compute TS value
289 295 double ts = 2.0 * (logL1 - m_logL0);
... ... @@ -431,7 +437,6 @@ void cttsmap::init_members(void)
431 437 m_publish = false;
432 438  
433 439 // Initialise protected members
434   - m_obs.clear();
435 440 m_binmin = -1;
436 441 m_binmax = -1;
437 442 m_logL0 = 0.0;
... ... @@ -463,11 +468,10 @@ void cttsmap::copy_members(const cttsmap&amp; app)
463 468 m_binmin = app.m_binmin;
464 469 m_binmax = app.m_binmax;
465 470 m_logL0 = app.m_logL0;
466   - m_obs = app.m_obs;
467 471 m_tsmap = app.m_tsmap;
  472 + m_statusmap = app.m_statusmap;
468 473 m_mapnames = app.m_mapnames;
469 474 m_maps = app.m_maps;
470   - m_statusmap = app.m_statusmap;
471 475  
472 476 // Clone protected members
473 477 m_testsource = (app.m_testsource != NULL) ? app.m_testsource->clone()
... ...
src/cttsmap/cttsmap.hpp
  View file @2aa78ba
... ... @@ -28,9 +28,7 @@
28 28 #define CTTSMAP_HPP
29 29  
30 30 /* __ Includes ___________________________________________________________ */
31   -#include "GammaLib.hpp"
32   -#include "GCTALib.hpp"
33   -#include "ctool.hpp"
  31 +#include "ctlikelihood.hpp"
34 32  
35 33 /* __Definitions _________________________________________________________ */
36 34 #define CTTSMAP_NAME "cttsmap"
... ... @@ -55,7 +53,7 @@
55 53 * grid of coordinates. The best fit results (TS, flux, index) are stored in
56 54 * maps which are saved in the output FITS files.
57 55 ***************************************************************************/
58   -class cttsmap : public ctool {
  56 +class cttsmap : public ctlikelihood {
59 57  
60 58 public:
61 59 // Constructors and destructors
... ... @@ -69,12 +67,11 @@ public:
69 67 cttsmap& operator=(const cttsmap& app);
70 68  
71 69 // Methods
72   - void clear(void);
73   - void run(void);
74   - void save(void);
75   - void publish(const std::string& name = "");
76   - const GObservations& obs(void) const;
77   - const GSkyMap& tsmap(void) const;
  70 + void clear(void);
  71 + void run(void);
  72 + void save(void);
  73 + void publish(const std::string& name = "");
  74 + const GSkyMap& tsmap(void) const;
78 75  
79 76 protected:
80 77 // Protected methods
... ... @@ -96,7 +93,6 @@ protected:
96 93 double m_logL0; //!< Likelihood value of null hypothesis
97 94  
98 95 // Protected members
99   - GObservations m_obs; //!< Observation container
100 96 GSkyMap m_tsmap; //!< TS map
101 97 GSkyMap m_statusmap; //!< Map of computed bins
102 98 std::vector<std::string> m_mapnames; //!< Names of free parameters
... ... @@ -106,18 +102,6 @@ protected:
106 102  
107 103  
108 104 /***********************************************************************//**
109   - * @brief Return observation container
110   - *
111   - * @return Reference to observation container
112   - ***************************************************************************/
113   -inline
114   -const GObservations& cttsmap::obs(void) const
115   -{
116   - return m_obs;
117   -}
118   -
119   -
120   -/***********************************************************************//**
121 105 * @brief Return TS skymap
122 106 *
123 107 * @return Reference to TS skymap
... ...
src/ctulimit/ctulimit.cpp
  View file @2aa78ba
... ... @@ -36,7 +36,6 @@
36 36 /* __ Method name definitions ____________________________________________ */
37 37 #define G_GET_MODEL_PARAMETER "ctulimit::get_model_parameter()"
38 38 #define G_UL_BISECTION "ctulimit::ul_bisection(double&, double&)"
39   -#define G_EVALUATE "ctulimit::evaluate(double&)"
40 39  
41 40 /* __ Debug definitions __________________________________________________ */
42 41  
... ... @@ -54,7 +53,7 @@
54 53 *
55 54 * Constructs an empty ctulimit tool.
56 55 ***************************************************************************/
57   -ctulimit::ctulimit(void) : ctool(CTULIMIT_NAME, CTULIMIT_VERSION)
  56 +ctulimit::ctulimit(void) : ctlikelihood(CTULIMIT_NAME, CTULIMIT_VERSION)
58 57 {
59 58 // Initialise members
60 59 init_members();
... ... @@ -69,18 +68,14 @@ ctulimit::ctulimit(void) : ctool(CTULIMIT_NAME, CTULIMIT_VERSION)
69 68 *
70 69 * param[in] obs Observation container.
71 70 *
72   - * Constructs an instance of the ctulimit tool that will operate on the
73   - * provided observation container.
  71 + * Constructs ctulimit tool from an observations container.
74 72 ***************************************************************************/
75 73 ctulimit::ctulimit(const GObservations& obs) :
76   - ctool(CTULIMIT_NAME, CTULIMIT_VERSION)
  74 + ctlikelihood(CTULIMIT_NAME, CTULIMIT_VERSION, obs)
77 75 {
78 76 // Initialise members
79 77 init_members();
80 78  
81   - // Set observations
82   - m_obs = obs;
83   -
84 79 // Return
85 80 return;
86 81 }
... ... @@ -97,7 +92,7 @@ ctulimit::ctulimit(const GObservations&amp; obs) :
97 92 * parameters that are provided as command line arguments.
98 93 ***************************************************************************/
99 94 ctulimit::ctulimit(int argc, char *argv[]) :
100   - ctool(CTULIMIT_NAME, CTULIMIT_VERSION, argc, argv)
  95 + ctlikelihood(CTULIMIT_NAME, CTULIMIT_VERSION, argc, argv)
101 96 {
102 97 // Initialise members
103 98 init_members();
... ... @@ -115,7 +110,7 @@ ctulimit::ctulimit(int argc, char *argv[]) :
115 110 * Constructs an instance of the ctulimit tool by copying information from
116 111 * another ctulimit tool.
117 112 ***************************************************************************/
118   -ctulimit::ctulimit(const ctulimit& app) : ctool(app)
  113 +ctulimit::ctulimit(const ctulimit& app) : ctlikelihood(app)
119 114 {
120 115 // Initialise members
121 116 init_members();
... ... @@ -163,7 +158,7 @@ ctulimit&amp; ctulimit::operator=(const ctulimit&amp; app)
163 158 if (this != &app) {
164 159  
165 160 // Copy base class members
166   - this->ctool::operator=(app);
  161 + this->ctlikelihood::operator=(app);
167 162  
168 163 // Free members
169 164 free_members();
... ... @@ -196,6 +191,7 @@ void ctulimit::clear(void)
196 191 {
197 192 // Free members
198 193 free_members();
  194 + this->ctlikelihood::free_members();
199 195 this->ctool::free_members();
200 196  
201 197 // Clear base class (needed to conserve tool name and version)
... ... @@ -203,6 +199,7 @@ void ctulimit::clear(void)
203 199  
204 200 // Initialise members
205 201 this->ctool::init_members();
  202 + this->ctlikelihood::init_members();
206 203 init_members();
207 204  
208 205 // Write header into logger
... ... @@ -378,8 +375,6 @@ void ctulimit::init_members(void)
378 375 m_chatter = static_cast<GChatter>(2);
379 376  
380 377 // Initialise protected members
381   - m_obs.clear();
382   - m_opt.clear();
383 378 m_dlogL = 0.0;
384 379 m_skymodel = NULL;
385 380 m_model_par = NULL;
... ... @@ -418,13 +413,11 @@ void ctulimit::copy_members(const ctulimit&amp; app)
418 413 m_chatter = app.m_chatter;
419 414  
420 415 // Copy protected members
421   - m_obs = app.m_obs;
422 416 m_dlogL = app.m_dlogL;
423 417 m_best_logL = app.m_best_logL;
424 418 m_diff_ulimit = app.m_diff_ulimit;
425 419 m_flux_ulimit = app.m_flux_ulimit;
426 420 m_eflux_ulimit = app.m_eflux_ulimit;
427   - m_opt = app.m_opt;
428 421  
429 422 // Extract model parameter
430 423 get_model_parameter();
... ... @@ -609,7 +602,7 @@ void ctulimit::ulimit_bisection(const double&amp; min, const double&amp; max)
609 602 double mid = (wrk_min + wrk_max) / 2.0;
610 603  
611 604 // Calculate function value
612   - double eval_mid = evaluate(mid);
  605 + double eval_mid = evaluate(*m_model_par, mid) - (m_best_logL + m_dlogL);
613 606  
614 607 // Check for convergence inside tolerance
615 608 if (std::abs(eval_mid) < m_tol) {
... ... @@ -633,53 +626,3 @@ void ctulimit::ulimit_bisection(const double&amp; min, const double&amp; max)
633 626 return;
634 627  
635 628 }
636   -
637   -
638   -/***********************************************************************//**
639   - * @brief Evaluates the log-likelihood
640   - *
641   - * @param[in] value Parameter factor value
642   - * @return Log-likelihood value
643   - *
644   - * Evaluates the log-likelihood as a function of the parameter of interest.
645   - ***************************************************************************/
646   -double ctulimit::evaluate(const double& value)
647   -{
648   - // Initialise log-likelihood value
649   - double logL = 0.0;
650   -
651   - // Check if given parameter is within boundaries
652   - if (value > m_model_par->factor_min() && value < m_model_par->factor_max()) {
653   -
654   - // Change parameter factor
655   - m_model_par->factor_value(value);
656   -
657   - // Fix parameter
658   - m_model_par->fix();
659   -
660   - // Re-optimize
661   - m_obs.optimize(m_opt);
662   -
663   - // Retrieve likelihood
664   - logL = m_obs.logL();
665   -
666   - // Free parameter
667   - m_model_par->free();
668   -
669   - } // endif: value was inside allowed range
670   -
671   - // ... otherwise signal that the parameter went outside the boundaries
672   - else {
673   - std::string msg = "Value of parameter \""+m_model_par->name()+"\" "
674   - "outside of validity range requested. To omit "
675   - "this error please enlarge the parameter range "
676   - "in the model XML file.";
677   - throw GException::invalid_value(G_EVALUATE, msg);
678   - }
679   -
680   - // Compute function value
681   - double logL_difference = logL - m_best_logL - m_dlogL;
682   -
683   - // Return
684   - return logL_difference;
685   -}
... ...
src/ctulimit/ctulimit.hpp
  View file @2aa78ba
... ... @@ -28,9 +28,7 @@
28 28 #define CTULIMIT_HPP
29 29  
30 30 /* __ Includes ___________________________________________________________ */
31   -#include "GammaLib.hpp"
32   -#include "GCTALib.hpp"
33   -#include "ctool.hpp"
  31 +#include "ctlikelihood.hpp"
34 32  
35 33 /* __Definitions _________________________________________________________ */
36 34 #define CTULIMIT_NAME "ctulimit"
... ... @@ -50,7 +48,7 @@
50 48 * During the computation the likelihood function is inspected to find the
51 49 * best value for the upper limit.
52 50 ***************************************************************************/
53   -class ctulimit : public ctool {
  51 +class ctulimit : public ctlikelihood {
54 52  
55 53 public:
56 54 // Constructors and destructors
... ... @@ -64,14 +62,12 @@ public:
64 62 ctulimit& operator=(const ctulimit& app);
65 63  
66 64 // Methods
67   - void clear(void);
68   - void run(void);
69   - void save(void);
70   - const GObservations& obs(void) const;
71   - const GOptimizer* opt(void) const;
72   - const double& diff_ulimit(void) const;
73   - const double& flux_ulimit(void) const;
74   - const double& eflux_ulimit(void) const;
  65 + void clear(void);
  66 + void run(void);
  67 + void save(void);
  68 + const double& diff_ulimit(void) const;
  69 + const double& flux_ulimit(void) const;
  70 + const double& eflux_ulimit(void) const;
75 71  
76 72 protected:
77 73 // Protected methods
... ... @@ -81,7 +77,6 @@ protected:
81 77 void get_parameters(void);
82 78 void get_model_parameter(void);
83 79 void ulimit_bisection(const double& parmin, const double& parmax);
84   - double evaluate(const double& value);
85 80  
86 81 // User parameters
87 82 std::string m_srcname; //!< Name of source which is moved around
... ... @@ -97,7 +92,6 @@ protected:
97 92 GChatter m_chatter; //!< Chattiness
98 93  
99 94 // Protected members
100   - GObservations m_obs; //!< Observation container
101 95 double m_dlogL; //!< Likelihood difference for upper limit computation
102 96 GModelSky* m_skymodel; //!< Pointer to sky model
103 97 GModelPar* m_model_par; //!< Pointer to model parameter
... ... @@ -105,35 +99,10 @@ protected:
105 99 double m_diff_ulimit; //!< Differential upper limit value
106 100 double m_flux_ulimit; //!< Flux upper limit value
107 101 double m_eflux_ulimit; //!< Energy flux upper limits
108   - GOptimizerLM m_opt; //!< Optimizer
109 102 };
110 103  
111 104  
112 105 /***********************************************************************//**
113   - * @brief Return observation container
114   - *
115   - * @return Reference to observation container.
116   - ***************************************************************************/
117   -inline
118   -const GObservations& ctulimit::obs(void) const
119   -{
120   - return m_obs;
121   -}
122   -
123   -
124   -/***********************************************************************//**
125   - * @brief Return optimizer
126   - *
127   - * @return Pointer to optimizer.
128   - ***************************************************************************/
129   -inline
130   -const GOptimizer* ctulimit::opt(void) const
131   -{
132   - return &m_opt;
133   -}
134   -
135   -
136   -/***********************************************************************//**
137 106 * @brief Return differential upper limit
138 107 *
139 108 * @return Differential upper flux limit.
... ...
src/support/Makefile.am
  View file @2aa78ba
... ... @@ -20,7 +20,7 @@
20 20 #############################################################################
21 21  
22 22 # Define sources for this directory
23   -sources = ctool.cpp ctool.hpp
  23 +sources = ctool.cpp ctool.hpp ctlikelihood.cpp ctlikelihood.hpp
24 24  
25 25 # Build libtool library
26 26 noinst_LTLIBRARIES = libsupport.la
... ...
src/support/ctlikelihood.cpp 0 → 100644
  View file @2aa78ba
  1 +/***************************************************************************
  2 + * ctlikelihood - Base class for likelihood tools *
  3 + * ----------------------------------------------------------------------- *
  4 + * copyright (C) 2016 by Juergen Knoedlseder *
  5 + * ----------------------------------------------------------------------- *
  6 + * *
  7 + * This program is free software: you can redistribute it and/or modify *
  8 + * it under the terms of the GNU General Public License as published by *
  9 + * the Free Software Foundation, either version 3 of the License, or *
  10 + * (at your option) any later version. *
  11 + * *
  12 + * This program is distributed in the hope that it will be useful, *
  13 + * but WITHOUT ANY WARRANTY; without even the implied warranty of *
  14 + * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the *
  15 + * GNU General Public License for more details. *
  16 + * *
  17 + * You should have received a copy of the GNU General Public License *
  18 + * along with this program. If not, see <http://www.gnu.org/licenses/>. *
  19 + * *
  20 + ***************************************************************************/
  21 +/**
  22 + * @file ctlikelihood.cpp
  23 + * @brief Likelihood tool base class implementation
  24 + * @author Juergen Knoedlseder
  25 + */
  26 +
  27 +/* __ Includes ___________________________________________________________ */
  28 +#ifdef HAVE_CONFIG_H
  29 +#include <config.h>
  30 +#endif
  31 +#include "GammaLib.hpp"
  32 +#include "GCTALib.hpp"
  33 +#include "ctlikelihood.hpp"
  34 +
  35 +//#include <iostream>
  36 +//#include "ctlikelihood.hpp"
  37 +//#include "GTools.hpp"
  38 +//#include "GOptimizer.hpp"
  39 +
  40 +/* __ Method name definitions ____________________________________________ */
  41 +#define G_EVALUATE "ctlikelihood::evaluate(GModelPar&, double&)"
  42 +
  43 +/* __ Debug definitions __________________________________________________ */
  44 +
  45 +/* __ Coding definitions _________________________________________________ */
  46 +
  47 +
  48 +/*==========================================================================
  49 + = =
  50 + = Constructors/destructors =
  51 + = =
  52 + ==========================================================================*/
  53 +
  54 +/***********************************************************************//**
  55 + * @brief Void constructor
  56 + *
  57 + * Constructs an empty likelihood tool.
  58 + ***************************************************************************/
  59 +ctlikelihood::ctlikelihood(void) : ctool()
  60 +{
  61 + // Initialise members
  62 + init_members();
  63 +
  64 + // Return
  65 + return;
  66 +}
  67 +
  68 +
  69 +/***********************************************************************//**
  70 + * @brief Name constructor
  71 + *
  72 + * @param[in] name Application name.
  73 + * @param[in] version Application version.
  74 + *
  75 + * Constructs a likelihood tool from the @p name and @p version. See the
  76 + * equivalent ctool constructor for details.
  77 + ***************************************************************************/
  78 +ctlikelihood::ctlikelihood(const std::string& name,
  79 + const std::string& version) : ctool(name, version)
  80 +{
  81 + // Initialise members
  82 + init_members();
  83 +
  84 + // Return
  85 + return;
  86 +}
  87 +
  88 +
  89 +/***********************************************************************//**
  90 + * @brief Observations constructor
  91 + *
  92 + * @param[in] name Application name.
  93 + * @param[in] version Application version.
  94 + * param[in] obs Observation container.
  95 + *
  96 + * Constructs a likelihood tool from the @p name, @p version and an
  97 + * observation container.
  98 + ***************************************************************************/
  99 +ctlikelihood::ctlikelihood(const std::string& name,
  100 + const std::string& version,
  101 + const GObservations& obs) : ctool(name, version)
  102 +{
  103 + // Initialise members
  104 + init_members();
  105 +
  106 + // Set observations
  107 + m_obs = obs;
  108 +
  109 + // Return
  110 + return;
  111 +}
  112 +
  113 +
  114 +/***********************************************************************//**
  115 + * @brief Command line constructor
  116 + *
  117 + * @param[in] name Application name.
  118 + * @param[in] version Application version.
  119 + * @param[in] argc Number of arguments in command line.
  120 + * @param[in] argv Array of command line arguments.
  121 + *
  122 + * Constructs a likelihood tool from the @p name, @p version and command
  123 + * line arguments. See the equivalent ctool constructor for details.
  124 + ***************************************************************************/
  125 +ctlikelihood::ctlikelihood(const std::string& name,
  126 + const std::string& version,
  127 + int argc,
  128 + char *argv[]) : ctool(name, version, argc, argv)
  129 +{
  130 + // Initialise members
  131 + init_members();
  132 +
  133 + // Return
  134 + return;
  135 +}
  136 +
  137 +
  138 +/***********************************************************************//**
  139 + * @brief Copy constructor
  140 + *
  141 + * @param[in] app Application.
  142 + *
  143 + * Constructs an instance of a likelihood tool by copying information from
  144 + * another likelihood tool.
  145 + ***************************************************************************/
  146 +ctlikelihood::ctlikelihood(const ctlikelihood& app) : ctool(app)
  147 +{
  148 + // Initialise members
  149 + init_members();
  150 +
  151 + // Copy members
  152 + copy_members(app);
  153 +
  154 + // Return
  155 + return;
  156 +}
  157 +
  158 +
  159 +/***********************************************************************//**
  160 + * @brief Destructor
  161 + *
  162 + * Destructs the ctlikelihood tool.
  163 + ***************************************************************************/
  164 +ctlikelihood::~ctlikelihood(void)
  165 +{
  166 + // Free members
  167 + free_members();
  168 +
  169 + // Return
  170 + return;
  171 +}
  172 +
  173 +
  174 +/*==========================================================================
  175 + = =
  176 + = Operators =
  177 + = =
  178 + ==========================================================================*/
  179 +
  180 +/***********************************************************************//**
  181 + * @brief Assignment operator
  182 + *
  183 + * @param[in] app Application.
  184 + * @return Application.
  185 + *
  186 + * Assigns a likelihood tool.
  187 + ***************************************************************************/
  188 +ctlikelihood& ctlikelihood::operator=(const ctlikelihood& app)
  189 +{
  190 + // Execute only if object is not identical
  191 + if (this != &app) {
  192 +
  193 + // Copy base class members
  194 + this->ctool::operator=(app);
  195 +
  196 + // Free members
  197 + free_members();
  198 +
  199 + // Initialise members
  200 + init_members();
  201 +
  202 + // Copy members
  203 + copy_members(app);
  204 +
  205 + } // endif: object was not identical
  206 +
  207 + // Return this object
  208 + return *this;
  209 +}
  210 +
  211 +
  212 +/*==========================================================================
  213 + = =
  214 + = Public methods =
  215 + = =
  216 + ==========================================================================*/
  217 +
  218 +
  219 +/*==========================================================================
  220 + = =
  221 + = Private methods =
  222 + = =
  223 + ==========================================================================*/
  224 +
  225 +/***********************************************************************//**
  226 + * @brief Initialise class members
  227 + ***************************************************************************/
  228 +void ctlikelihood::init_members(void)
  229 +{
  230 + // Initialise members
  231 + m_obs.clear();
  232 + m_opt.clear();
  233 +
  234 + // Return
  235 + return;
  236 +}
  237 +
  238 +
  239 +/***********************************************************************//**
  240 + * @brief Copy class members
  241 + *
  242 + * @param[in] app Application.
  243 + ***************************************************************************/
  244 +void ctlikelihood::copy_members(const ctlikelihood& app)
  245 +{
  246 + // Copy members
  247 + m_obs = app.m_obs;
  248 + m_opt = app.m_opt;
  249 +
  250 + // Return
  251 + return;
  252 +}
  253 +
  254 +
  255 +/***********************************************************************//**
  256 + * @brief Delete class members
  257 + ***************************************************************************/
  258 +void ctlikelihood::free_members(void)
  259 +{
  260 + // Return
  261 + return;
  262 +}
  263 +
  264 +
  265 +/***********************************************************************//**
  266 + * @brief Evaluates the log-likelihood function
  267 + *
  268 + * @param[in] par Model parameter
  269 + * @param[in] value Model parameter factor value
  270 + * @return Log-likelihood function
  271 + *
  272 + * Evaluates the log-likelihood function at a given @p value.
  273 + ***************************************************************************/
  274 +double ctlikelihood::evaluate(GModelPar& par, const double& value)
  275 +{
  276 + // Initialise log-likelihood value
  277 + double logL = 0.0;
  278 +
  279 + // Throw an exception if the parameter is below the minimum boundary
  280 + if (value < par.factor_min()) {
  281 + std::string msg = "Value "+gammalib::str(value)+" of parameter \""+
  282 + par.name()+"\" is below its minimum boundary "+
  283 + gammalib::str(par.factor_min())+". To omit this "
  284 + "error please lower the minimum parameter boundary.";
  285 + throw GException::invalid_value(G_EVALUATE, msg);
  286 + }
  287 +
  288 + // Throw an exception if the parameter is above the maximum boundary
  289 + if (value > par.factor_max()) {
  290 + std::string msg = "Value "+gammalib::str(value)+" of parameter \""+
  291 + par.name()+"\" is above its maximum boundary "+
  292 + gammalib::str(par.factor_max())+". To omit this "
  293 + "error please raise the maximum parameter boundary.";
  294 + throw GException::invalid_value(G_EVALUATE, msg);
  295 + }
  296 +
  297 + // Change parameter factor
  298 + par.factor_value(value);
  299 +
  300 + // Fix parameter
  301 + par.fix();
  302 +
  303 + // Re-optimize log-likelihood
  304 + m_obs.optimize(m_opt);
  305 +
  306 + // Free parameter
  307 + par.free();
  308 +
  309 + // Retrieve log-likelihood
  310 + logL = m_obs.logL();
  311 +
  312 + // Return log-likelihood
  313 + return logL;
  314 +}
... ...
src/support/ctlikelihood.hpp 0 → 100644
  View file @2aa78ba
  1 +/***************************************************************************
  2 + * ctlikelihood - Base class for likelihood tools *
  3 + * ----------------------------------------------------------------------- *
  4 + * copyright (C) 2016 by Juergen Knoedlseder *
  5 + * ----------------------------------------------------------------------- *
  6 + * *
  7 + * This program is free software: you can redistribute it and/or modify *
  8 + * it under the terms of the GNU General Public License as published by *
  9 + * the Free Software Foundation, either version 3 of the License, or *
  10 + * (at your option) any later version. *
  11 + * *
  12 + * This program is distributed in the hope that it will be useful, *
  13 + * but WITHOUT ANY WARRANTY; without even the implied warranty of *
  14 + * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the *
  15 + * GNU General Public License for more details. *
  16 + * *
  17 + * You should have received a copy of the GNU General Public License *
  18 + * along with this program. If not, see <http://www.gnu.org/licenses/>. *
  19 + * *
  20 + ***************************************************************************/
  21 +/**
  22 + * @file ctlikelihood.hpp
  23 + * @brief Likelihood tool base class interface definition
  24 + * @author Juergen Knoedlseder
  25 + */
  26 +
  27 +#ifndef CTLIKELIHOOD_HPP
  28 +#define CTLIKELIHOOD_HPP
  29 +
  30 +/* __ Includes ___________________________________________________________ */
  31 +#include "ctool.hpp"
  32 +
  33 +/* __Definitions _________________________________________________________ */
  34 +
  35 +
  36 +/***********************************************************************//**
  37 + * @class ctlikelihood
  38 + *
  39 + * @brief Base class for likelihood tools
  40 + ***************************************************************************/
  41 +class ctlikelihood : public ctool {
  42 +
  43 +public:
  44 + // Constructors and destructors
  45 + ctlikelihood(void);
  46 + ctlikelihood(const std::string& name, const std::string& version);
  47 + ctlikelihood(const std::string& name, const std::string& version,
  48 + const GObservations& obs);
  49 + ctlikelihood(const std::string& name, const std::string& version,
  50 + int argc, char* argv[]);
  51 + ctlikelihood(const ctlikelihood& app);
  52 + virtual ~ctlikelihood(void);
  53 +
  54 + // Operators
  55 + ctlikelihood& operator=(const ctlikelihood& app);
  56 +
  57 + // Pure virtual methods
  58 + virtual void clear(void) = 0;
  59 + virtual void run(void) = 0;
  60 + virtual void save(void) = 0;
  61 +
  62 + // Methods
  63 + const GObservations& obs(void) const;
  64 + const GOptimizer* opt(void) const;
  65 +
  66 +protected:
  67 + // Protected methods
  68 + void init_members(void);
  69 + void copy_members(const ctlikelihood& app);
  70 + void free_members(void);
  71 + double evaluate(GModelPar& par, const double& value);
  72 +
  73 + // Protected members
  74 + GObservations m_obs; //!< Observation container
  75 + GOptimizerLM m_opt; //!< Optimizer
  76 +};
  77 +
  78 +
  79 +/***********************************************************************//**
  80 + * @brief Return observation container
  81 + *
  82 + * @return Reference to observation container.
  83 + *
  84 + * Returns a reference to the observation container.
  85 + ***************************************************************************/
  86 +inline
  87 +const GObservations& ctlikelihood::obs(void) const
  88 +{
  89 + return m_obs;
  90 +}
  91 +
  92 +
  93 +/***********************************************************************//**
  94 + * @brief Return optimizer
  95 + *
  96 + * @return Pointer to optimizer.
  97 + *
  98 + * Returns a pointer to the optimizer.
  99 + ***************************************************************************/
  100 +inline
  101 +const GOptimizer* ctlikelihood::opt(void) const
  102 +{
  103 + return &m_opt;
  104 +}
  105 +
  106 +#endif /* CTLIKELIHOOD_HPP */
... ...
src/support/ctool.cpp
  View file @2aa78ba
... ... @@ -1655,8 +1655,8 @@ GEnergies ctool::insert_energy_boundaries(const GEnergies&amp; energies,
1655 1655 /***********************************************************************//**
1656 1656 * @brief Determine the counts cube layer usage.
1657 1657 *
1658   - * @param[in] ebounds_cube Energy boundaries of the counts cube.
1659   - * @param[in] ebounds_list Energy boundaries of the event list.
  1658 + * @param[in] cube_ebounds Energy boundaries of the counts cube.
  1659 + * @param[in] list_ebounds Energy boundaries of the event list.
1660 1660 * @return Vector of usage flags.
1661 1661 *
1662 1662 * Determines a vector of counts cube layer usage flags that signal whether
... ...
src/support/ctool.hpp
  View file @2aa78ba
... ... @@ -37,7 +37,13 @@
37 37 /***********************************************************************//**
38 38 * @class ctool
39 39 *
40   - * @brief ctool base class
  40 + * @brief Base class for ctools
  41 + *
  42 + * This class is the generic base class for all ctools. It enforces
  43 + * implementation of the clear(), run() and save() methods by all ctools,
  44 + * and it implements the execute() method.
  45 + *
  46 + * The class also provides a number of helper methods for ctools.
41 47 ***************************************************************************/
42 48 class ctool : public GApplication {
43 49 public:
... ...
test/test_ctulimit.py
  View file @2aa78ba
... ... @@ -103,9 +103,9 @@ class Test(test):
103 103 ulimit = ctools.ctulimit()
104 104  
105 105 # Check that empty ctulimit tool holds zero upper limits
106   - self.test_value(ulimit.diff_ulimit(), 0.0, 1.0e-23,
  106 + self.test_value(ulimit.diff_ulimit(), 0.0, 1.0e-21,
107 107 'Check differential upper limit')
108   - self.test_value(ulimit.flux_ulimit(), 0.0, 1.0e-17,
  108 + self.test_value(ulimit.flux_ulimit(), 0.0, 1.0e-16,
109 109 'Check upper limit on photon flux')
110 110 self.test_value(ulimit.eflux_ulimit(), 0.0, 1.0e-16,
111 111 'Check upper limit on energy flux')
... ... @@ -132,9 +132,9 @@ class Test(test):
132 132 ulimit.save()
133 133  
134 134 # Check results
135   - self.test_value(ulimit.diff_ulimit(), 8.86803e-18, 1.0e-23,
  135 + self.test_value(ulimit.diff_ulimit(), 8.86803e-18, 1.0e-21,
136 136 'Check differential upper limit')
137   - self.test_value(ulimit.flux_ulimit(), 6.10509e-12, 1.0e-17,
  137 + self.test_value(ulimit.flux_ulimit(), 6.10509e-12, 1.0e-16,
138 138 'Check upper limit on photon flux')
139 139 self.test_value(ulimit.eflux_ulimit(), 2.75669e-11, 1.0e-16,
140 140 'Check upper limit on energy flux')
... ... @@ -150,10 +150,10 @@ class Test(test):
150 150 cpy_ulimit = ulimit.copy()
151 151  
152 152 # Check results of copy
153   - self.test_value(cpy_ulimit.diff_ulimit(), 8.86803e-18, 1.0e-23,
  153 + self.test_value(cpy_ulimit.diff_ulimit(), 8.86803e-18, 1.0e-21,
154 154 'Check differential upper limit')
155   - self.test_value(cpy_ulimit.flux_ulimit(), 6.10509e-12, 1.0e-17,
156   - 'ulimit upper limit on photon flux')
  155 + self.test_value(cpy_ulimit.flux_ulimit(), 6.10509e-12, 1.0e-16,
  156 + 'Check upper limit on photon flux')
157 157 self.test_value(cpy_ulimit.eflux_ulimit(), 2.75669e-11, 1.0e-16,
158 158 'Check upper limit on energy flux')
159 159  
... ... @@ -165,10 +165,10 @@ class Test(test):
165 165 cpy_ulimit.execute()
166 166  
167 167 # Check results
168   - self.test_value(cpy_ulimit.diff_ulimit(), 8.86803e-18, 1.0e-23,
  168 + self.test_value(cpy_ulimit.diff_ulimit(), 8.86803e-18, 1.0e-21,
169 169 'Check differential upper limit')
170   - self.test_value(cpy_ulimit.flux_ulimit(), 6.10509e-12, 1.0e-17,
171   - 'ulimit upper limit on photon flux')
  170 + self.test_value(cpy_ulimit.flux_ulimit(), 6.10509e-12, 1.0e-16,
  171 + 'Check upper limit on photon flux')
172 172 self.test_value(cpy_ulimit.eflux_ulimit(), 2.75669e-11, 1.0e-16,
173 173 'Check upper limit on energy flux')
174 174  
... ... @@ -179,9 +179,9 @@ class Test(test):
179 179 cpy_ulimit.clear()
180 180  
181 181 # Check that empty ctulimit tool holds zero upper limits
182   - self.test_value(cpy_ulimit.diff_ulimit(), 0.0, 1.0e-23,
  182 + self.test_value(cpy_ulimit.diff_ulimit(), 0.0, 1.0e-21,
183 183 'Check differential upper limit')
184   - self.test_value(cpy_ulimit.flux_ulimit(), 0.0, 1.0e-17,
  184 + self.test_value(cpy_ulimit.flux_ulimit(), 0.0, 1.0e-16,
185 185 'Check upper limit on photon flux')
186 186 self.test_value(cpy_ulimit.eflux_ulimit(), 0.0, 1.0e-16,
187 187 'Check upper limit on energy flux')
... ... @@ -190,11 +190,7 @@ class Test(test):
190 190 # time. This should avoid the necessity to recompute the maximum
191 191 # likelihood
192 192 ulimit = ctools.ctulimit(obs)
193   - #ulimit['inobs'] = self._events
194   - #ulimit['inmodel'] = self._model
195 193 ulimit['srcname'] = 'Crab'
196   - #ulimit['caldb'] = self._caldb
197   - #ulimit['irf'] = self._irf
198 194 ulimit['logfile'] = 'ctulimit_py3.log'
199 195 ulimit['chatter'] = 4
200 196  
... ... @@ -203,10 +199,10 @@ class Test(test):
203 199 ulimit.execute()
204 200  
205 201 # Check results
206   - self.test_value(ulimit.diff_ulimit(), 8.86803e-18, 1.0e-23,
  202 + self.test_value(ulimit.diff_ulimit(), 8.86803e-18, 1.0e-21,
207 203 'Check differential upper limit')
208   - self.test_value(ulimit.flux_ulimit(), 6.10509e-12, 1.0e-17,
209   - 'ulimit upper limit on photon flux')
  204 + self.test_value(ulimit.flux_ulimit(), 6.10509e-12, 1.0e-16,
  205 + 'Check upper limit on photon flux')
210 206 self.test_value(ulimit.eflux_ulimit(), 2.75669e-11, 1.0e-16,
211 207 'Check upper limit on energy flux')
212 208  
... ...