diff --git a/ChangeLog b/ChangeLog
index b8867aa..edfbd95 100644
--- a/ChangeLog
+++ b/ChangeLog
@@ -3,6 +3,7 @@
* Version 1.2.0 released
========================
+ Introduce ctlikelihood base class (#1828)
Remove ctselect::time_reference() method
Add ctool and cscript pars() methods (#1838)
Fix ctbin::save() bug for empty observation container (#1835)
diff --git a/NEWS b/NEWS
index 3540a36..e4b3436 100644
--- a/NEWS
+++ b/NEWS
@@ -24,7 +24,10 @@ all parameters in a dictionary (#1838).
Support model
-------------
-None
+The ctlikelihood base class was introduced from which all likelihood ctools
+now derive (ctlike, cttsmap, cterror, ctulimit). This avoids code duplication
+and provides all tools with a standard interface. All tools now have the
+obs() and opt() method (#1828).
ctobssim - CTA observation simulation
@@ -57,12 +60,16 @@ ctlike - Maximum likelihood fitting
The clear() method now calls GApplication::clear() to conserve the properties
of the tool (#1835).
+The class now derives from ctlikelihood (#1828).
+
cttsmap - Generation of Test Statistic map
------------------------------------------
The clear() method now calls GApplication::clear() to conserve the properties
of the tool (#1835).
+The class now derives from ctlikelihood (#1828).
+
ctskymap - CTA sky mapping
--------------------------
@@ -134,12 +141,16 @@ of the tool (#1835).
The tool now also supports spectral models with the new naming conventions, and
in particular the "PhotonFlux" and "EnergyFlux" names (#1832).
+The class now derives from ctlikelihood (#1828).
+
cterror - Likelihood profile error computation
----------------------------------------------
The clear() method now calls GApplication::clear() to conserve the properties
of the tool (#1835).
+The class now derives from ctlikelihood (#1828).
+
cscripts
--------
diff --git a/pyext/Makefile.am b/pyext/Makefile.am
index 770b777..69dd74f 100755
--- a/pyext/Makefile.am
+++ b/pyext/Makefile.am
@@ -168,8 +168,8 @@ CLEANFILES = $(CLEAN_FILES) $(DEPFILES)
# Files to be removed in build directory
clean-local:
rm -rf build
- rm -rf $(builddir)/ctools/*.py
- rm -rf $(builddir)/ctools/*.pyc
+ rm -rf $(builddir)/ctools/c*.py
+ rm -rf $(builddir)/ctools/c*.pyc
rm -rf $(builddir)/ctools/*_wrap.cpp
rm -rf $(builddir)/ctools/_tools.so
rm -rf $(builddir)/ctools/__pycache__
diff --git a/pyext/cterror.i b/pyext/cterror.i
index 5ac43b2..89dd73e 100644
--- a/pyext/cterror.i
+++ b/pyext/cterror.i
@@ -1,7 +1,7 @@
/***************************************************************************
* cterror - Parameter error calculation tool *
* ----------------------------------------------------------------------- *
- * copyright (C) 2015 by Florent Forest *
+ * copyright (C) 2015-2016 by Florent Forest *
* ----------------------------------------------------------------------- *
* *
* This program is free software: you can redistribute it and/or modify *
@@ -34,9 +34,8 @@
* @class cterror
*
* @brief Parameter error calculation tool
- *
***************************************************************************/
-class cterror : public ctool {
+class cterror : public ctlikelihood {
public:
// Constructors and destructors
@@ -47,11 +46,9 @@ public:
virtual ~cterror(void);
// Methods
- void clear(void);
- void run(void);
- void save(void);
- const GObservations& obs(void) const;
- const GOptimizer* opt(void) const;
+ void clear(void);
+ void run(void);
+ void save(void);
};
diff --git a/pyext/ctlike.i b/pyext/ctlike.i
index def8e83..7fb4d52 100644
--- a/pyext/ctlike.i
+++ b/pyext/ctlike.i
@@ -1,7 +1,7 @@
/***************************************************************************
* ctlike - Maximum likelihood fitting tool *
* ----------------------------------------------------------------------- *
- * copyright (C) 2010-2014 by Jurgen Knodlseder *
+ * copyright (C) 2010-2016 by Jurgen Knodlseder *
* ----------------------------------------------------------------------- *
* *
* This program is free software: you can redistribute it and/or modify *
@@ -34,7 +34,7 @@
*
* @brief Maximum likelihood fitting tool
***************************************************************************/
-class ctlike : public ctool {
+class ctlike : public ctlikelihood {
public:
// Constructors and destructors
ctlike(void);
@@ -44,11 +44,9 @@ public:
virtual ~ctlike(void);
// Methods
- void clear(void);
- void run(void);
- void save(void);
- const GObservations& obs(void) const;
- const GOptimizer* opt(void) const;
+ void clear(void);
+ void run(void);
+ void save(void);
};
diff --git a/pyext/ctlikelihood.i b/pyext/ctlikelihood.i
new file mode 100644
index 0000000..47450ae
--- /dev/null
+++ b/pyext/ctlikelihood.i
@@ -0,0 +1,66 @@
+/***************************************************************************
+ * ctlikelihood - Base class for likelihood tools *
+ * ----------------------------------------------------------------------- *
+ * copyright (C) 2016 by Juergen Knoedlseder *
+ * ----------------------------------------------------------------------- *
+ * *
+ * This program is free software: you can redistribute it and/or modify *
+ * it under the terms of the GNU General Public License as published by *
+ * the Free Software Foundation, either version 3 of the License, or *
+ * (at your option) any later version. *
+ * *
+ * This program is distributed in the hope that it will be useful, *
+ * but WITHOUT ANY WARRANTY; without even the implied warranty of *
+ * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the *
+ * GNU General Public License for more details. *
+ * *
+ * You should have received a copy of the GNU General Public License *
+ * along with this program. If not, see . *
+ * *
+ ***************************************************************************/
+/**
+ * @file ctlikelihood.i
+ * @brief Likelihood tool base class interface definition
+ * @author Juergen Knoedlseder
+ */
+
+%{
+/* Put headers and other declarations here that are needed for compilation */
+#include "ctlikelihood.hpp"
+%}
+
+
+/***********************************************************************//**
+ * @class ctlikelihood
+ *
+ * @brief Base class for likelihood tools
+ ***************************************************************************/
+class ctlikelihood : public ctool {
+
+public:
+ // Constructors and destructors
+ ctlikelihood(void);
+ ctlikelihood(const std::string& name, const std::string& version);
+ ctlikelihood(const std::string& name, const std::string& version,
+ const GObservations& obs);
+ ctlikelihood(const std::string& name, const std::string& version,
+ int argc, char* argv[]);
+ ctlikelihood(const ctlikelihood& app);
+ virtual ~ctlikelihood(void);
+
+ // Pure virtual methods
+ virtual void clear(void) = 0;
+ virtual void run(void) = 0;
+ virtual void save(void) = 0;
+
+ // Methods
+ const GObservations& obs(void) const;
+ const GOptimizer* opt(void) const;
+};
+
+
+/***********************************************************************//**
+ * @brief Likelihood tool Python extensions
+ ***************************************************************************/
+%extend ctlikelihood {
+}
diff --git a/pyext/ctool.i b/pyext/ctool.i
index 05e5324..5e85e0d 100644
--- a/pyext/ctool.i
+++ b/pyext/ctool.i
@@ -62,7 +62,7 @@ public:
/***********************************************************************//**
* @class ctool
*
- * @brief ctool base class
+ * @brief Base class for ctools
***************************************************************************/
class ctool : public GApplication {
public:
diff --git a/pyext/ctools/tools.i b/pyext/ctools/tools.i
index b726f62..1b3fc82 100755
--- a/pyext/ctools/tools.i
+++ b/pyext/ctools/tools.i
@@ -44,8 +44,11 @@
%import(module="gammalib.base") "GBase.i";
%import(module="gammalib.app") "GApplication.i";
-/* __ ctools _____________________________________________________________ */
+/* __ ctools base classes ________________________________________________ */
%include "ctool.i"
+%include "ctlikelihood.i"
+
+/* __ ctools _____________________________________________________________ */
%include "ctbin.i"
%include "ctobssim.i"
%include "ctlike.i"
diff --git a/pyext/cttsmap.i b/pyext/cttsmap.i
index 0954d2c..dd64342 100644
--- a/pyext/cttsmap.i
+++ b/pyext/cttsmap.i
@@ -34,7 +34,8 @@
*
* @brief TS map calculation tool
***************************************************************************/
-class cttsmap : public ctool {
+class cttsmap : public ctlikelihood {
+
public:
// Constructors and destructors
cttsmap(void);
@@ -44,12 +45,11 @@ public:
virtual ~cttsmap(void);
// Methods
- void clear(void);
- void run(void);
- void save(void);
- void publish(const std::string& name = "");
- const GObservations& obs(void) const;
- const GSkyMap& tsmap(void) const;
+ void clear(void);
+ void run(void);
+ void save(void);
+ void publish(const std::string& name = "");
+ const GSkyMap& tsmap(void) const;
};
diff --git a/pyext/ctulimit.i b/pyext/ctulimit.i
index 289189d..135b27b 100644
--- a/pyext/ctulimit.i
+++ b/pyext/ctulimit.i
@@ -1,7 +1,7 @@
/***************************************************************************
* ctulimit - Upper limit calculation tool *
* ----------------------------------------------------------------------- *
- * copyright (C) 2015 by Michael Mayer *
+ * copyright (C) 2015-2016 by Michael Mayer *
* ----------------------------------------------------------------------- *
* *
* This program is free software: you can redistribute it and/or modify *
@@ -34,7 +34,7 @@
*
* @brief Upper limit calculation tool
***************************************************************************/
-class ctulimit : public ctool {
+class ctulimit : public ctlikelihood {
public:
// Constructors and destructors
ctulimit(void);
@@ -44,14 +44,12 @@ public:
virtual ~ctulimit(void);
// Methods
- void clear(void);
- void run(void);
- void save(void);
- const GObservations& obs(void) const;
- const GOptimizer* opt(void) const;
- const double& diff_ulimit(void) const;
- const double& flux_ulimit(void) const;
- const double& eflux_ulimit(void) const;
+ void clear(void);
+ void run(void);
+ void save(void);
+ const double& diff_ulimit(void) const;
+ const double& flux_ulimit(void) const;
+ const double& eflux_ulimit(void) const;
};
diff --git a/src/cterror/cterror.cpp b/src/cterror/cterror.cpp
index 0ef808b..d7b1dad 100644
--- a/src/cterror/cterror.cpp
+++ b/src/cterror/cterror.cpp
@@ -35,7 +35,6 @@
/* __ Method name definitions ____________________________________________ */
#define G_ERR_BISECTION "cterror::error_bisection(double&, double&)"
-#define G_EVALUATE "cterror::evaluate(double&)"
/* __ Debug definitions __________________________________________________ */
@@ -53,7 +52,7 @@
*
* Constructs an empty cterror tool.
***************************************************************************/
-cterror::cterror(void) : ctool(CTERROR_NAME, CTERROR_VERSION)
+cterror::cterror(void) : ctlikelihood(CTERROR_NAME, CTERROR_VERSION)
{
// Initialise members
init_members();
@@ -68,24 +67,19 @@ cterror::cterror(void) : ctool(CTERROR_NAME, CTERROR_VERSION)
*
* param[in] obs Observation container.
*
- * Constructs an instance of the cterror tool that will operate on the
- * provided observation container.
+ * Constructs cterror tool from an observations container.
***************************************************************************/
cterror::cterror(const GObservations& obs) :
- ctool(CTERROR_NAME, CTERROR_VERSION)
+ ctlikelihood(CTERROR_NAME, CTERROR_VERSION, obs)
{
// Initialise members
init_members();
- // Set observations
- m_obs = obs;
-
// Return
return;
}
-
/***********************************************************************//**
* @brief Command line constructor
*
@@ -96,7 +90,7 @@ cterror::cterror(const GObservations& obs) :
* parameters that are provided as command line arguments.
***************************************************************************/
cterror::cterror(int argc, char *argv[]) :
- ctool(CTERROR_NAME, CTERROR_VERSION, argc, argv)
+ ctlikelihood(CTERROR_NAME, CTERROR_VERSION, argc, argv)
{
// Initialise members
init_members();
@@ -114,7 +108,7 @@ cterror::cterror(int argc, char *argv[]) :
* Constructs an instance of the cterror tool by copying information from
* another ctulimit tool.
***************************************************************************/
-cterror::cterror(const cterror& app) : ctool(app)
+cterror::cterror(const cterror& app) : ctlikelihood(app)
{
// Initialise members
init_members();
@@ -162,7 +156,7 @@ cterror& cterror::operator=(const cterror& app)
if (this != &app) {
// Copy base class members
- this->ctool::operator=(app);
+ this->ctlikelihood::operator=(app);
// Free members
free_members();
@@ -195,6 +189,7 @@ void cterror::clear(void)
{
// Free members
free_members();
+ this->ctlikelihood::free_members();
this->ctool::free_members();
// Clear base class (needed to conserve tool name and version)
@@ -202,6 +197,7 @@ void cterror::clear(void)
// Initialise members
this->ctool::init_members();
+ this->ctlikelihood::init_members();
init_members();
// Write header into logger
@@ -384,19 +380,17 @@ void cterror::init_members(void)
// Initialise user parameters
m_srcname.clear();
m_outmodel.clear();
- m_confidence = 0.68;
- m_tol = 1.0e-3;
- m_max_iter = 50;
- m_value = 0.0;
- m_chatter = static_cast(2);
+ m_confidence = 0.68;
+ m_tol = 1.0e-3;
+ m_max_iter = 50;
+ m_apply_edisp = false;
+ m_chatter = static_cast(2);
// Initialise protected members
- m_obs.clear();
- m_opt.clear();
+ m_value = 0.0;
m_dlogL = 0.0;
m_best_logL = 0.0;
m_model_par = NULL;
- m_apply_edisp = false;
// Set optimizer parameters
m_opt.max_iter(m_max_iter);
@@ -424,10 +418,9 @@ void cterror::copy_members(const cterror& app)
m_chatter = app.m_chatter;
// Copy protected members
- m_obs = app.m_obs;
+ m_value = app.m_value;
m_dlogL = app.m_dlogL;
m_best_logL = app.m_best_logL;
- m_opt = app.m_opt;
m_model_par = NULL;
// Return
@@ -548,7 +541,7 @@ double cterror::error_bisection(const double& min, const double& max)
mid = (wrk_min + wrk_max) / 2.0;
// Calculate function value
- double eval_mid = evaluate(mid);
+ double eval_mid = evaluate(*m_model_par, mid) - (m_best_logL + m_dlogL);
// Write interval into logger
log_value(EXPLICIT, " Iteration "+gammalib::str(iter),
@@ -597,53 +590,3 @@ double cterror::error_bisection(const double& min, const double& max)
return mid;
}
-
-
-/***********************************************************************//**
- * @brief Evaluates the log-likelihood
- *
- * @param[in] value Parameter factor value
- * @return Log-likelihood value
- *
- * Evaluates the log-likelihood as a function of the parameter of interest.
- ***************************************************************************/
-double cterror::evaluate(const double& value)
-{
- // Initialise log-likelihood value
- double logL = 0.0;
-
- // Check if given parameter is within boundaries
- if (value > m_model_par->factor_min() && value < m_model_par->factor_max()) {
-
- // Change parameter factor
- m_model_par->factor_value(value);
-
- // Fix parameter
- m_model_par->fix();
-
- // Re-optimize
- m_obs.optimize(m_opt);
-
- // Retrieve likelihood
- logL = m_obs.logL();
-
- // Free parameter
- m_model_par->free();
-
- } // endif: value was inside allowed range
-
- // ... otherwise signal that the parameter went outside the boundaries
- else {
- std::string msg = "Value of parameter \""+m_model_par->name()+"\" "
- "outside of validity range requested. To omit "
- "this error please enlarge the parameter range "
- "in the model XML file.";
- throw GException::invalid_value(G_EVALUATE, msg);
- }
-
- // Compute function value
- double logL_difference = logL - m_best_logL - m_dlogL;
-
- // Return
- return logL_difference;
-}
diff --git a/src/cterror/cterror.hpp b/src/cterror/cterror.hpp
index 249666b..24293f0 100644
--- a/src/cterror/cterror.hpp
+++ b/src/cterror/cterror.hpp
@@ -28,10 +28,7 @@
#define CTERROR_HPP
/* __ Includes ___________________________________________________________ */
-#include
-#include "GammaLib.hpp"
-#include "GCTALib.hpp"
-#include "ctool.hpp"
+#include "ctlikelihood.hpp"
/* __Definitions _________________________________________________________ */
#define CTERROR_NAME "cterror"
@@ -43,7 +40,7 @@
*
* @brief Parameter error calculation tool
***************************************************************************/
-class cterror : public ctool {
+class cterror : public ctlikelihood {
public:
// Constructors and destructors
@@ -57,11 +54,9 @@ public:
cterror& operator=(const cterror& app);
// Methods
- void clear(void);
- void run(void);
- void save(void);
- const GObservations& obs(void) const;
- const GOptimizer* opt(void) const;
+ void clear(void);
+ void run(void);
+ void save(void);
protected:
// Protected methods
@@ -69,7 +64,6 @@ protected:
void copy_members(const cterror& app);
void free_members(void);
void get_parameters(void);
- double evaluate(const double& value);
double error_bisection(const double& min, const double& max);
// User parameters
@@ -78,40 +72,14 @@ protected:
double m_confidence; //!< Confidence level
double m_tol; //!< Tolerance for limit determination
int m_max_iter; //!< Maximum number of iterations
- double m_value; //!< Parameter value
bool m_apply_edisp; //!< Apply energy dispersion?
GChatter m_chatter; //!< Chattiness
// Protected members
- GObservations m_obs; //!< Observation container
+ double m_value; //!< Parameter value
double m_dlogL; //!< Likelihood difference for upper limit computation
GModelPar* m_model_par; //!< Pointer to model parameter
double m_best_logL; //!< Best fit log likelihood of given model
- GOptimizerLM m_opt; //!< Optimizer
};
-
-/***********************************************************************//**
- * @brief Return observation container
- *
- * @return Reference to observation container.
- ***************************************************************************/
-inline
-const GObservations& cterror::obs(void) const
-{
- return m_obs;
-}
-
-
-/***********************************************************************//**
- * @brief Return optimizer
- *
- * @return Pointer to optimizer.
- ***************************************************************************/
-inline
-const GOptimizer* cterror::opt(void) const
-{
- return &m_opt;
-}
-
#endif /* CTERROR_HPP */
diff --git a/src/ctlike/ctlike.cpp b/src/ctlike/ctlike.cpp
index ec566be..e03a43a 100644
--- a/src/ctlike/ctlike.cpp
+++ b/src/ctlike/ctlike.cpp
@@ -48,7 +48,7 @@
/***********************************************************************//**
* @brief Void constructor
***************************************************************************/
-ctlike::ctlike(void) : ctool(CTLIKE_NAME, CTLIKE_VERSION)
+ctlike::ctlike(void) : ctlikelihood(CTLIKE_NAME, CTLIKE_VERSION)
{
// Initialise members
init_members();
@@ -63,17 +63,14 @@ ctlike::ctlike(void) : ctool(CTLIKE_NAME, CTLIKE_VERSION)
*
* param[in] obs Observation container.
*
- * This method creates an instance of the class by copying an existing
- * observations container.
+ * Constructs ctlike tool from an observations container.
***************************************************************************/
-ctlike::ctlike(const GObservations& obs) : ctool(CTLIKE_NAME, CTLIKE_VERSION)
+ctlike::ctlike(const GObservations& obs) :
+ ctlikelihood(CTLIKE_NAME, CTLIKE_VERSION, obs)
{
// Initialise members
init_members();
- // Set observations
- m_obs = obs;
-
// Return
return;
}
@@ -84,9 +81,12 @@ ctlike::ctlike(const GObservations& obs) : ctool(CTLIKE_NAME, CTLIKE_VERSION)
*
* @param[in] argc Number of arguments in command line.
* @param[in] argv Array of command line arguments.
+ *
+ * Constructs an instance of the ctlike tool that will parse user parameters
+ * that are provided as command line arguments.
***************************************************************************/
ctlike::ctlike(int argc, char *argv[]) :
- ctool(CTLIKE_NAME, CTLIKE_VERSION, argc, argv)
+ ctlikelihood(CTLIKE_NAME, CTLIKE_VERSION, argc, argv)
{
// Initialise members
init_members();
@@ -101,7 +101,7 @@ ctlike::ctlike(int argc, char *argv[]) :
*
* @param[in] app Application.
***************************************************************************/
-ctlike::ctlike(const ctlike& app) : ctool(app)
+ctlike::ctlike(const ctlike& app) : ctlikelihood(app)
{
// Initialise members
init_members();
@@ -145,7 +145,7 @@ ctlike& ctlike::operator=(const ctlike& app)
if (this != &app) {
// Copy base class members
- this->ctool::operator=(app);
+ this->ctlikelihood::operator=(app);
// Free members
free_members();
@@ -178,6 +178,7 @@ void ctlike::clear(void)
{
// Free members
free_members();
+ this->ctlikelihood::free_members();
this->ctool::free_members();
// Clear base class (needed to conserve tool name and version)
@@ -185,6 +186,7 @@ void ctlike::clear(void)
// Initialise members
this->ctool::init_members();
+ this->ctlikelihood::init_members();
init_members();
// Write header into logger
@@ -291,7 +293,7 @@ void ctlike::run(void)
// Write results into logger
log_header1(NORMAL, "Maximum likelihood optimisation results");
- log_string(NORMAL, m_opt->print(m_chatter));
+ log_string(NORMAL, m_opt.print(m_chatter));
log_value(NORMAL, "Maximum log likelihood", gammalib::str(m_logL,3));
log_value(NORMAL, "Observed events (Nobs)", gammalib::str(num_events,3));
log_value(NORMAL, "Predicted events (Npred)", gammalib::str(npred,3)+
@@ -398,10 +400,10 @@ void ctlike::get_parameters(void)
// Set optimizer logger
if (logNormal()) {
- static_cast(m_opt)->logger(&log);
+ static_cast(&m_opt)->logger(&log);
}
else {
- static_cast(m_opt)->logger(NULL);
+ static_cast(&m_opt)->logger(NULL);
}
// Write parameters into logger
@@ -442,7 +444,7 @@ void ctlike::optimize_lm(void)
}
// Perform LM optimization
- m_obs.optimize(*m_opt);
+ m_obs.optimize(m_opt);
// Optionally refit
if (m_refit) {
@@ -453,7 +455,7 @@ void ctlike::optimize_lm(void)
log.indent(1);
// Optimise again
- m_obs.optimize(*m_opt);
+ m_obs.optimize(m_opt);
}
@@ -464,10 +466,10 @@ void ctlike::optimize_lm(void)
(m_obs.function()).curvature())->print());
// Compute errors
- m_obs.errors(*m_opt);
+ m_obs.errors(m_opt);
// Store maximum log likelihood value
- m_logL = -(m_opt->value());
+ m_logL = -(m_opt.value());
// Remove indent
log.indent(0);
@@ -490,7 +492,7 @@ double ctlike::reoptimize_lm(void)
log.indent(1);
// Create a clone of the optimizer for the re-optimisation
- GOptimizer* opt = m_opt->clone();
+ GOptimizer* opt = m_opt.clone();
// Perform LM optimization
m_obs.optimize(*opt);
@@ -526,12 +528,10 @@ void ctlike::init_members(void)
{
// Initialise members
m_outmodel.clear();
- m_obs.clear();
m_refit = false;
m_max_iter = 100; // Set maximum number of iterations
m_max_stall = 10; // Set maximum number of stalls
m_logL = 0.0;
- m_opt = NULL;
m_apply_edisp = false;
m_fix_spat_for_ts = false;
m_chatter = static_cast(2);
@@ -539,15 +539,9 @@ void ctlike::init_members(void)
// Set logger properties
log.date(true);
- // Allocate LM optimizer
- GOptimizerLM* opt = new GOptimizerLM();
-
// Set optimizer parameters
- opt->max_iter(m_max_iter);
- opt->max_stalls(m_max_stall);
-
- // Set optimizer pointer
- m_opt = opt;
+ m_opt.max_iter(m_max_iter);
+ m_opt.max_stalls(m_max_stall);
// Return
return;
@@ -564,11 +558,9 @@ void ctlike::copy_members(const ctlike& app)
// Copy attributes
m_refit = app.m_refit;
m_outmodel = app.m_outmodel;
- m_obs = app.m_obs;
m_max_iter = app.m_max_iter;
m_max_stall = app.m_max_stall;
m_logL = app.m_logL;
- m_opt = app.m_opt->clone();
m_apply_edisp = app.m_apply_edisp;
m_fix_spat_for_ts = app.m_fix_spat_for_ts;
m_chatter = app.m_chatter;
@@ -583,12 +575,6 @@ void ctlike::copy_members(const ctlike& app)
***************************************************************************/
void ctlike::free_members(void)
{
- // Free members
- if (m_opt != NULL) delete m_opt;
-
- // Mark pointers as free
- m_opt = NULL;
-
// Return
return;
}
diff --git a/src/ctlike/ctlike.hpp b/src/ctlike/ctlike.hpp
index a88767a..11e99cc 100644
--- a/src/ctlike/ctlike.hpp
+++ b/src/ctlike/ctlike.hpp
@@ -28,14 +28,11 @@
#define CTLIKE_HPP
/* __ Includes ___________________________________________________________ */
-#include
-#include "GammaLib.hpp"
-#include "GCTALib.hpp"
-#include "ctool.hpp"
+#include "ctlikelihood.hpp"
/* __Definitions _________________________________________________________ */
#define CTLIKE_NAME "ctlike"
-#define CTLIKE_VERSION "1.1.0"
+#define CTLIKE_VERSION "1.2.0"
/***********************************************************************//**
@@ -43,7 +40,7 @@
*
* @brief Maximum likelihood fitting tool
***************************************************************************/
-class ctlike : public ctool {
+class ctlike : public ctlikelihood {
public:
// Constructors and destructors
@@ -57,11 +54,9 @@ public:
ctlike& operator=(const ctlike& app);
// Methods
- void clear(void);
- void run(void);
- void save(void);
- const GObservations& obs(void) const;
- const GOptimizer* opt(void) const;
+ void clear(void);
+ void run(void);
+ void save(void);
protected:
// Protected methods
@@ -80,35 +75,9 @@ protected:
GChatter m_chatter; //!< Chattiness
// Members
- GObservations m_obs; //!< Observations
int m_max_iter; //!< Maximum number of iterations
int m_max_stall; //!< Maximum number of stalls
double m_logL; //!< Maximum log likelihood
- GOptimizer* m_opt; //!< Optimizer
};
-
-/***********************************************************************//**
- * @brief Return observation container
- *
- * @return Reference to observation container
- ***************************************************************************/
-inline
-const GObservations& ctlike::obs(void) const
-{
- return m_obs;
-}
-
-
-/***********************************************************************//**
- * @brief Return optimizer
- *
- * @return Pointer to optimizer
- ***************************************************************************/
-inline
-const GOptimizer* ctlike::opt(void) const
-{
- return m_opt;
-}
-
#endif /* CTLIKE_HPP */
diff --git a/src/ctmodel/ctmodel.cpp b/src/ctmodel/ctmodel.cpp
index 7c1d6b2..c3b6c0b 100644
--- a/src/ctmodel/ctmodel.cpp
+++ b/src/ctmodel/ctmodel.cpp
@@ -598,7 +598,7 @@ void ctmodel::get_parameters(void)
* "expcube", "psfcube", and "bkgcube" parameters). In that case, it appended
* a dummy event cube to the observation.
*
- * @TODO Support stacked energy dispersion
+ * @todo Support stacked energy dispersion
***************************************************************************/
void ctmodel::get_obs(void)
{
diff --git a/src/ctobssim/ctobssim.cpp b/src/ctobssim/ctobssim.cpp
index 66136bd..aafd088 100644
--- a/src/ctobssim/ctobssim.cpp
+++ b/src/ctobssim/ctobssim.cpp
@@ -1039,10 +1039,11 @@ void ctobssim::simulate_interval(GCTAObservation* obs,
* @brief Simulate source events for a time slice
*
* @param[in] obs Pointer on CTA observation.
+ * @param[in] rsp Pointer on CTA response.
* @param[in,out] events Pointer on CTA event list.
- * @param[in] models Model list.
- * @param[in] tmin Start time.
- * @param[in] tmax Stop time.
+ * @param[in] model Model.
+ * @param[in] tstart Start time.
+ * @param[in] tstop Stop time.
* @param[in] etrue_min Minimum true energy.
* @param[in] etrue_max Maximum true energy.
* @param[in] ereco_min Minimum reconstructed energy.
diff --git a/src/cttsmap/cttsmap.cpp b/src/cttsmap/cttsmap.cpp
index 664d983..b4e6165 100644
--- a/src/cttsmap/cttsmap.cpp
+++ b/src/cttsmap/cttsmap.cpp
@@ -48,8 +48,10 @@
/***********************************************************************//**
* @brief Void constructor
+ *
+ * Constructs an empty cttsmap tool.
***************************************************************************/
-cttsmap::cttsmap(void) : ctool(CTTSMAP_NAME, CTTSMAP_VERSION)
+cttsmap::cttsmap(void) : ctlikelihood(CTTSMAP_NAME, CTTSMAP_VERSION)
{
// Initialise members
init_members();
@@ -64,32 +66,30 @@ cttsmap::cttsmap(void) : ctool(CTTSMAP_NAME, CTTSMAP_VERSION)
*
* param[in] obs Observation container.
*
- * This method creates an instance of the class by copying an existing
- * observations container.
+ * Constructs cttsmap tool from an observations container.
***************************************************************************/
cttsmap::cttsmap(const GObservations& obs) :
- ctool(CTTSMAP_NAME, CTTSMAP_VERSION)
+ ctlikelihood(CTTSMAP_NAME, CTTSMAP_VERSION, obs)
{
// Initialise members
init_members();
- // Set observations
- m_obs = obs;
-
// Return
return;
}
-
/***********************************************************************//**
* @brief Command line constructor
*
* @param[in] argc Number of arguments in command line.
* @param[in] argv Array of command line arguments.
+ *
+ * Constructs an instance of the cttsmap tool that will parse user parameters
+ * that are provided as command line arguments.
***************************************************************************/
cttsmap::cttsmap(int argc, char *argv[]) :
- ctool(CTTSMAP_NAME, CTTSMAP_VERSION, argc, argv)
+ ctlikelihood(CTTSMAP_NAME, CTTSMAP_VERSION, argc, argv)
{
// Initialise members
init_members();
@@ -104,7 +104,7 @@ cttsmap::cttsmap(int argc, char *argv[]) :
*
* @param[in] app Application.
***************************************************************************/
-cttsmap::cttsmap(const cttsmap& app) : ctool(app)
+cttsmap::cttsmap(const cttsmap& app) : ctlikelihood(app)
{
// Initialise members
init_members();
@@ -148,7 +148,7 @@ cttsmap& cttsmap::operator=(const cttsmap& app)
if (this != &app) {
// Copy base class members
- this->ctool::operator=(app);
+ this->ctlikelihood::operator=(app);
// Free members
free_members();
@@ -181,6 +181,7 @@ void cttsmap::clear(void)
{
// Free members
free_members();
+ this->ctlikelihood::free_members();
this->ctool::free_members();
// Clear base class (needed to conserve tool name and version)
@@ -188,6 +189,7 @@ void cttsmap::clear(void)
// Initialise members
this->ctool::init_members();
+ this->ctlikelihood::init_members();
init_members();
// Write header into logger
@@ -232,9 +234,13 @@ void cttsmap::run(void)
int binmin = (m_binmin == -1) ? 0 : m_binmin;
int binmax = (m_binmax == -1) ? m_tsmap.npix() : m_binmax;
- // Initialise optimizer
- GOptimizerLM* opt = (logExplicit()) ? new GOptimizerLM(&log)
- : new GOptimizerLM();
+ // Set optimizer logger
+ if (logExplicit()) {
+ static_cast(&m_opt)->logger(&log);
+ }
+ else {
+ static_cast(&m_opt)->logger(NULL);
+ }
// Store initial models
GModels models_orig = m_obs.models();
@@ -253,8 +259,8 @@ void cttsmap::run(void)
// Compute likelihood without the test source
m_obs.models(models);
- m_obs.optimize(*opt);
- m_logL0 = -(opt->value());
+ m_obs.optimize(m_opt);
+ m_logL0 = -(m_opt.value());
}
@@ -280,10 +286,10 @@ void cttsmap::run(void)
m_obs.models(models);
// Optimize observation container
- m_obs.optimize(*opt);
+ m_obs.optimize(m_opt);
// Retrieve the Likelihood value
- double logL1 = -(opt->value());
+ double logL1 = -(m_opt.value());
// Compute TS value
double ts = 2.0 * (logL1 - m_logL0);
@@ -431,7 +437,6 @@ void cttsmap::init_members(void)
m_publish = false;
// Initialise protected members
- m_obs.clear();
m_binmin = -1;
m_binmax = -1;
m_logL0 = 0.0;
@@ -463,11 +468,10 @@ void cttsmap::copy_members(const cttsmap& app)
m_binmin = app.m_binmin;
m_binmax = app.m_binmax;
m_logL0 = app.m_logL0;
- m_obs = app.m_obs;
m_tsmap = app.m_tsmap;
+ m_statusmap = app.m_statusmap;
m_mapnames = app.m_mapnames;
m_maps = app.m_maps;
- m_statusmap = app.m_statusmap;
// Clone protected members
m_testsource = (app.m_testsource != NULL) ? app.m_testsource->clone()
diff --git a/src/cttsmap/cttsmap.hpp b/src/cttsmap/cttsmap.hpp
index 9a21805..7f2e12c 100644
--- a/src/cttsmap/cttsmap.hpp
+++ b/src/cttsmap/cttsmap.hpp
@@ -28,9 +28,7 @@
#define CTTSMAP_HPP
/* __ Includes ___________________________________________________________ */
-#include "GammaLib.hpp"
-#include "GCTALib.hpp"
-#include "ctool.hpp"
+#include "ctlikelihood.hpp"
/* __Definitions _________________________________________________________ */
#define CTTSMAP_NAME "cttsmap"
@@ -55,7 +53,7 @@
* grid of coordinates. The best fit results (TS, flux, index) are stored in
* maps which are saved in the output FITS files.
***************************************************************************/
-class cttsmap : public ctool {
+class cttsmap : public ctlikelihood {
public:
// Constructors and destructors
@@ -69,12 +67,11 @@ public:
cttsmap& operator=(const cttsmap& app);
// Methods
- void clear(void);
- void run(void);
- void save(void);
- void publish(const std::string& name = "");
- const GObservations& obs(void) const;
- const GSkyMap& tsmap(void) const;
+ void clear(void);
+ void run(void);
+ void save(void);
+ void publish(const std::string& name = "");
+ const GSkyMap& tsmap(void) const;
protected:
// Protected methods
@@ -96,7 +93,6 @@ protected:
double m_logL0; //!< Likelihood value of null hypothesis
// Protected members
- GObservations m_obs; //!< Observation container
GSkyMap m_tsmap; //!< TS map
GSkyMap m_statusmap; //!< Map of computed bins
std::vector m_mapnames; //!< Names of free parameters
@@ -106,18 +102,6 @@ protected:
/***********************************************************************//**
- * @brief Return observation container
- *
- * @return Reference to observation container
- ***************************************************************************/
-inline
-const GObservations& cttsmap::obs(void) const
-{
- return m_obs;
-}
-
-
-/***********************************************************************//**
* @brief Return TS skymap
*
* @return Reference to TS skymap
diff --git a/src/ctulimit/ctulimit.cpp b/src/ctulimit/ctulimit.cpp
index 654e69b..0cde83e 100644
--- a/src/ctulimit/ctulimit.cpp
+++ b/src/ctulimit/ctulimit.cpp
@@ -36,7 +36,6 @@
/* __ Method name definitions ____________________________________________ */
#define G_GET_MODEL_PARAMETER "ctulimit::get_model_parameter()"
#define G_UL_BISECTION "ctulimit::ul_bisection(double&, double&)"
-#define G_EVALUATE "ctulimit::evaluate(double&)"
/* __ Debug definitions __________________________________________________ */
@@ -54,7 +53,7 @@
*
* Constructs an empty ctulimit tool.
***************************************************************************/
-ctulimit::ctulimit(void) : ctool(CTULIMIT_NAME, CTULIMIT_VERSION)
+ctulimit::ctulimit(void) : ctlikelihood(CTULIMIT_NAME, CTULIMIT_VERSION)
{
// Initialise members
init_members();
@@ -69,18 +68,14 @@ ctulimit::ctulimit(void) : ctool(CTULIMIT_NAME, CTULIMIT_VERSION)
*
* param[in] obs Observation container.
*
- * Constructs an instance of the ctulimit tool that will operate on the
- * provided observation container.
+ * Constructs ctulimit tool from an observations container.
***************************************************************************/
ctulimit::ctulimit(const GObservations& obs) :
- ctool(CTULIMIT_NAME, CTULIMIT_VERSION)
+ ctlikelihood(CTULIMIT_NAME, CTULIMIT_VERSION, obs)
{
// Initialise members
init_members();
- // Set observations
- m_obs = obs;
-
// Return
return;
}
@@ -97,7 +92,7 @@ ctulimit::ctulimit(const GObservations& obs) :
* parameters that are provided as command line arguments.
***************************************************************************/
ctulimit::ctulimit(int argc, char *argv[]) :
- ctool(CTULIMIT_NAME, CTULIMIT_VERSION, argc, argv)
+ ctlikelihood(CTULIMIT_NAME, CTULIMIT_VERSION, argc, argv)
{
// Initialise members
init_members();
@@ -115,7 +110,7 @@ ctulimit::ctulimit(int argc, char *argv[]) :
* Constructs an instance of the ctulimit tool by copying information from
* another ctulimit tool.
***************************************************************************/
-ctulimit::ctulimit(const ctulimit& app) : ctool(app)
+ctulimit::ctulimit(const ctulimit& app) : ctlikelihood(app)
{
// Initialise members
init_members();
@@ -163,7 +158,7 @@ ctulimit& ctulimit::operator=(const ctulimit& app)
if (this != &app) {
// Copy base class members
- this->ctool::operator=(app);
+ this->ctlikelihood::operator=(app);
// Free members
free_members();
@@ -196,6 +191,7 @@ void ctulimit::clear(void)
{
// Free members
free_members();
+ this->ctlikelihood::free_members();
this->ctool::free_members();
// Clear base class (needed to conserve tool name and version)
@@ -203,6 +199,7 @@ void ctulimit::clear(void)
// Initialise members
this->ctool::init_members();
+ this->ctlikelihood::init_members();
init_members();
// Write header into logger
@@ -378,8 +375,6 @@ void ctulimit::init_members(void)
m_chatter = static_cast(2);
// Initialise protected members
- m_obs.clear();
- m_opt.clear();
m_dlogL = 0.0;
m_skymodel = NULL;
m_model_par = NULL;
@@ -418,13 +413,11 @@ void ctulimit::copy_members(const ctulimit& app)
m_chatter = app.m_chatter;
// Copy protected members
- m_obs = app.m_obs;
m_dlogL = app.m_dlogL;
m_best_logL = app.m_best_logL;
m_diff_ulimit = app.m_diff_ulimit;
m_flux_ulimit = app.m_flux_ulimit;
m_eflux_ulimit = app.m_eflux_ulimit;
- m_opt = app.m_opt;
// Extract model parameter
get_model_parameter();
@@ -609,7 +602,7 @@ void ctulimit::ulimit_bisection(const double& min, const double& max)
double mid = (wrk_min + wrk_max) / 2.0;
// Calculate function value
- double eval_mid = evaluate(mid);
+ double eval_mid = evaluate(*m_model_par, mid) - (m_best_logL + m_dlogL);
// Check for convergence inside tolerance
if (std::abs(eval_mid) < m_tol) {
@@ -633,53 +626,3 @@ void ctulimit::ulimit_bisection(const double& min, const double& max)
return;
}
-
-
-/***********************************************************************//**
- * @brief Evaluates the log-likelihood
- *
- * @param[in] value Parameter factor value
- * @return Log-likelihood value
- *
- * Evaluates the log-likelihood as a function of the parameter of interest.
- ***************************************************************************/
-double ctulimit::evaluate(const double& value)
-{
- // Initialise log-likelihood value
- double logL = 0.0;
-
- // Check if given parameter is within boundaries
- if (value > m_model_par->factor_min() && value < m_model_par->factor_max()) {
-
- // Change parameter factor
- m_model_par->factor_value(value);
-
- // Fix parameter
- m_model_par->fix();
-
- // Re-optimize
- m_obs.optimize(m_opt);
-
- // Retrieve likelihood
- logL = m_obs.logL();
-
- // Free parameter
- m_model_par->free();
-
- } // endif: value was inside allowed range
-
- // ... otherwise signal that the parameter went outside the boundaries
- else {
- std::string msg = "Value of parameter \""+m_model_par->name()+"\" "
- "outside of validity range requested. To omit "
- "this error please enlarge the parameter range "
- "in the model XML file.";
- throw GException::invalid_value(G_EVALUATE, msg);
- }
-
- // Compute function value
- double logL_difference = logL - m_best_logL - m_dlogL;
-
- // Return
- return logL_difference;
-}
diff --git a/src/ctulimit/ctulimit.hpp b/src/ctulimit/ctulimit.hpp
index b202202..cb16e5f 100644
--- a/src/ctulimit/ctulimit.hpp
+++ b/src/ctulimit/ctulimit.hpp
@@ -28,9 +28,7 @@
#define CTULIMIT_HPP
/* __ Includes ___________________________________________________________ */
-#include "GammaLib.hpp"
-#include "GCTALib.hpp"
-#include "ctool.hpp"
+#include "ctlikelihood.hpp"
/* __Definitions _________________________________________________________ */
#define CTULIMIT_NAME "ctulimit"
@@ -50,7 +48,7 @@
* During the computation the likelihood function is inspected to find the
* best value for the upper limit.
***************************************************************************/
-class ctulimit : public ctool {
+class ctulimit : public ctlikelihood {
public:
// Constructors and destructors
@@ -64,14 +62,12 @@ public:
ctulimit& operator=(const ctulimit& app);
// Methods
- void clear(void);
- void run(void);
- void save(void);
- const GObservations& obs(void) const;
- const GOptimizer* opt(void) const;
- const double& diff_ulimit(void) const;
- const double& flux_ulimit(void) const;
- const double& eflux_ulimit(void) const;
+ void clear(void);
+ void run(void);
+ void save(void);
+ const double& diff_ulimit(void) const;
+ const double& flux_ulimit(void) const;
+ const double& eflux_ulimit(void) const;
protected:
// Protected methods
@@ -81,7 +77,6 @@ protected:
void get_parameters(void);
void get_model_parameter(void);
void ulimit_bisection(const double& parmin, const double& parmax);
- double evaluate(const double& value);
// User parameters
std::string m_srcname; //!< Name of source which is moved around
@@ -97,7 +92,6 @@ protected:
GChatter m_chatter; //!< Chattiness
// Protected members
- GObservations m_obs; //!< Observation container
double m_dlogL; //!< Likelihood difference for upper limit computation
GModelSky* m_skymodel; //!< Pointer to sky model
GModelPar* m_model_par; //!< Pointer to model parameter
@@ -105,35 +99,10 @@ protected:
double m_diff_ulimit; //!< Differential upper limit value
double m_flux_ulimit; //!< Flux upper limit value
double m_eflux_ulimit; //!< Energy flux upper limits
- GOptimizerLM m_opt; //!< Optimizer
};
/***********************************************************************//**
- * @brief Return observation container
- *
- * @return Reference to observation container.
- ***************************************************************************/
-inline
-const GObservations& ctulimit::obs(void) const
-{
- return m_obs;
-}
-
-
-/***********************************************************************//**
- * @brief Return optimizer
- *
- * @return Pointer to optimizer.
- ***************************************************************************/
-inline
-const GOptimizer* ctulimit::opt(void) const
-{
- return &m_opt;
-}
-
-
-/***********************************************************************//**
* @brief Return differential upper limit
*
* @return Differential upper flux limit.
diff --git a/src/support/Makefile.am b/src/support/Makefile.am
index 98f1333..4c00567 100644
--- a/src/support/Makefile.am
+++ b/src/support/Makefile.am
@@ -20,7 +20,7 @@
#############################################################################
# Define sources for this directory
-sources = ctool.cpp ctool.hpp
+sources = ctool.cpp ctool.hpp ctlikelihood.cpp ctlikelihood.hpp
# Build libtool library
noinst_LTLIBRARIES = libsupport.la
diff --git a/src/support/ctlikelihood.cpp b/src/support/ctlikelihood.cpp
new file mode 100644
index 0000000..78341d3
--- /dev/null
+++ b/src/support/ctlikelihood.cpp
@@ -0,0 +1,314 @@
+/***************************************************************************
+ * ctlikelihood - Base class for likelihood tools *
+ * ----------------------------------------------------------------------- *
+ * copyright (C) 2016 by Juergen Knoedlseder *
+ * ----------------------------------------------------------------------- *
+ * *
+ * This program is free software: you can redistribute it and/or modify *
+ * it under the terms of the GNU General Public License as published by *
+ * the Free Software Foundation, either version 3 of the License, or *
+ * (at your option) any later version. *
+ * *
+ * This program is distributed in the hope that it will be useful, *
+ * but WITHOUT ANY WARRANTY; without even the implied warranty of *
+ * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the *
+ * GNU General Public License for more details. *
+ * *
+ * You should have received a copy of the GNU General Public License *
+ * along with this program. If not, see . *
+ * *
+ ***************************************************************************/
+/**
+ * @file ctlikelihood.cpp
+ * @brief Likelihood tool base class implementation
+ * @author Juergen Knoedlseder
+ */
+
+/* __ Includes ___________________________________________________________ */
+#ifdef HAVE_CONFIG_H
+#include
+#endif
+#include "GammaLib.hpp"
+#include "GCTALib.hpp"
+#include "ctlikelihood.hpp"
+
+//#include
+//#include "ctlikelihood.hpp"
+//#include "GTools.hpp"
+//#include "GOptimizer.hpp"
+
+/* __ Method name definitions ____________________________________________ */
+#define G_EVALUATE "ctlikelihood::evaluate(GModelPar&, double&)"
+
+/* __ Debug definitions __________________________________________________ */
+
+/* __ Coding definitions _________________________________________________ */
+
+
+/*==========================================================================
+ = =
+ = Constructors/destructors =
+ = =
+ ==========================================================================*/
+
+/***********************************************************************//**
+ * @brief Void constructor
+ *
+ * Constructs an empty likelihood tool.
+ ***************************************************************************/
+ctlikelihood::ctlikelihood(void) : ctool()
+{
+ // Initialise members
+ init_members();
+
+ // Return
+ return;
+}
+
+
+/***********************************************************************//**
+ * @brief Name constructor
+ *
+ * @param[in] name Application name.
+ * @param[in] version Application version.
+ *
+ * Constructs a likelihood tool from the @p name and @p version. See the
+ * equivalent ctool constructor for details.
+ ***************************************************************************/
+ctlikelihood::ctlikelihood(const std::string& name,
+ const std::string& version) : ctool(name, version)
+{
+ // Initialise members
+ init_members();
+
+ // Return
+ return;
+}
+
+
+/***********************************************************************//**
+ * @brief Observations constructor
+ *
+ * @param[in] name Application name.
+ * @param[in] version Application version.
+ * param[in] obs Observation container.
+ *
+ * Constructs a likelihood tool from the @p name, @p version and an
+ * observation container.
+ ***************************************************************************/
+ctlikelihood::ctlikelihood(const std::string& name,
+ const std::string& version,
+ const GObservations& obs) : ctool(name, version)
+{
+ // Initialise members
+ init_members();
+
+ // Set observations
+ m_obs = obs;
+
+ // Return
+ return;
+}
+
+
+/***********************************************************************//**
+ * @brief Command line constructor
+ *
+ * @param[in] name Application name.
+ * @param[in] version Application version.
+ * @param[in] argc Number of arguments in command line.
+ * @param[in] argv Array of command line arguments.
+ *
+ * Constructs a likelihood tool from the @p name, @p version and command
+ * line arguments. See the equivalent ctool constructor for details.
+ ***************************************************************************/
+ctlikelihood::ctlikelihood(const std::string& name,
+ const std::string& version,
+ int argc,
+ char *argv[]) : ctool(name, version, argc, argv)
+{
+ // Initialise members
+ init_members();
+
+ // Return
+ return;
+}
+
+
+/***********************************************************************//**
+ * @brief Copy constructor
+ *
+ * @param[in] app Application.
+ *
+ * Constructs an instance of a likelihood tool by copying information from
+ * another likelihood tool.
+ ***************************************************************************/
+ctlikelihood::ctlikelihood(const ctlikelihood& app) : ctool(app)
+{
+ // Initialise members
+ init_members();
+
+ // Copy members
+ copy_members(app);
+
+ // Return
+ return;
+}
+
+
+/***********************************************************************//**
+ * @brief Destructor
+ *
+ * Destructs the ctlikelihood tool.
+ ***************************************************************************/
+ctlikelihood::~ctlikelihood(void)
+{
+ // Free members
+ free_members();
+
+ // Return
+ return;
+}
+
+
+/*==========================================================================
+ = =
+ = Operators =
+ = =
+ ==========================================================================*/
+
+/***********************************************************************//**
+ * @brief Assignment operator
+ *
+ * @param[in] app Application.
+ * @return Application.
+ *
+ * Assigns a likelihood tool.
+ ***************************************************************************/
+ctlikelihood& ctlikelihood::operator=(const ctlikelihood& app)
+{
+ // Execute only if object is not identical
+ if (this != &app) {
+
+ // Copy base class members
+ this->ctool::operator=(app);
+
+ // Free members
+ free_members();
+
+ // Initialise members
+ init_members();
+
+ // Copy members
+ copy_members(app);
+
+ } // endif: object was not identical
+
+ // Return this object
+ return *this;
+}
+
+
+/*==========================================================================
+ = =
+ = Public methods =
+ = =
+ ==========================================================================*/
+
+
+/*==========================================================================
+ = =
+ = Private methods =
+ = =
+ ==========================================================================*/
+
+/***********************************************************************//**
+ * @brief Initialise class members
+ ***************************************************************************/
+void ctlikelihood::init_members(void)
+{
+ // Initialise members
+ m_obs.clear();
+ m_opt.clear();
+
+ // Return
+ return;
+}
+
+
+/***********************************************************************//**
+ * @brief Copy class members
+ *
+ * @param[in] app Application.
+ ***************************************************************************/
+void ctlikelihood::copy_members(const ctlikelihood& app)
+{
+ // Copy members
+ m_obs = app.m_obs;
+ m_opt = app.m_opt;
+
+ // Return
+ return;
+}
+
+
+/***********************************************************************//**
+ * @brief Delete class members
+ ***************************************************************************/
+void ctlikelihood::free_members(void)
+{
+ // Return
+ return;
+}
+
+
+/***********************************************************************//**
+ * @brief Evaluates the log-likelihood function
+ *
+ * @param[in] par Model parameter
+ * @param[in] value Model parameter factor value
+ * @return Log-likelihood function
+ *
+ * Evaluates the log-likelihood function at a given @p value.
+ ***************************************************************************/
+double ctlikelihood::evaluate(GModelPar& par, const double& value)
+{
+ // Initialise log-likelihood value
+ double logL = 0.0;
+
+ // Throw an exception if the parameter is below the minimum boundary
+ if (value < par.factor_min()) {
+ std::string msg = "Value "+gammalib::str(value)+" of parameter \""+
+ par.name()+"\" is below its minimum boundary "+
+ gammalib::str(par.factor_min())+". To omit this "
+ "error please lower the minimum parameter boundary.";
+ throw GException::invalid_value(G_EVALUATE, msg);
+ }
+
+ // Throw an exception if the parameter is above the maximum boundary
+ if (value > par.factor_max()) {
+ std::string msg = "Value "+gammalib::str(value)+" of parameter \""+
+ par.name()+"\" is above its maximum boundary "+
+ gammalib::str(par.factor_max())+". To omit this "
+ "error please raise the maximum parameter boundary.";
+ throw GException::invalid_value(G_EVALUATE, msg);
+ }
+
+ // Change parameter factor
+ par.factor_value(value);
+
+ // Fix parameter
+ par.fix();
+
+ // Re-optimize log-likelihood
+ m_obs.optimize(m_opt);
+
+ // Free parameter
+ par.free();
+
+ // Retrieve log-likelihood
+ logL = m_obs.logL();
+
+ // Return log-likelihood
+ return logL;
+}
diff --git a/src/support/ctlikelihood.hpp b/src/support/ctlikelihood.hpp
new file mode 100644
index 0000000..1e600ae
--- /dev/null
+++ b/src/support/ctlikelihood.hpp
@@ -0,0 +1,106 @@
+/***************************************************************************
+ * ctlikelihood - Base class for likelihood tools *
+ * ----------------------------------------------------------------------- *
+ * copyright (C) 2016 by Juergen Knoedlseder *
+ * ----------------------------------------------------------------------- *
+ * *
+ * This program is free software: you can redistribute it and/or modify *
+ * it under the terms of the GNU General Public License as published by *
+ * the Free Software Foundation, either version 3 of the License, or *
+ * (at your option) any later version. *
+ * *
+ * This program is distributed in the hope that it will be useful, *
+ * but WITHOUT ANY WARRANTY; without even the implied warranty of *
+ * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the *
+ * GNU General Public License for more details. *
+ * *
+ * You should have received a copy of the GNU General Public License *
+ * along with this program. If not, see . *
+ * *
+ ***************************************************************************/
+/**
+ * @file ctlikelihood.hpp
+ * @brief Likelihood tool base class interface definition
+ * @author Juergen Knoedlseder
+ */
+
+#ifndef CTLIKELIHOOD_HPP
+#define CTLIKELIHOOD_HPP
+
+/* __ Includes ___________________________________________________________ */
+#include "ctool.hpp"
+
+/* __Definitions _________________________________________________________ */
+
+
+/***********************************************************************//**
+ * @class ctlikelihood
+ *
+ * @brief Base class for likelihood tools
+ ***************************************************************************/
+class ctlikelihood : public ctool {
+
+public:
+ // Constructors and destructors
+ ctlikelihood(void);
+ ctlikelihood(const std::string& name, const std::string& version);
+ ctlikelihood(const std::string& name, const std::string& version,
+ const GObservations& obs);
+ ctlikelihood(const std::string& name, const std::string& version,
+ int argc, char* argv[]);
+ ctlikelihood(const ctlikelihood& app);
+ virtual ~ctlikelihood(void);
+
+ // Operators
+ ctlikelihood& operator=(const ctlikelihood& app);
+
+ // Pure virtual methods
+ virtual void clear(void) = 0;
+ virtual void run(void) = 0;
+ virtual void save(void) = 0;
+
+ // Methods
+ const GObservations& obs(void) const;
+ const GOptimizer* opt(void) const;
+
+protected:
+ // Protected methods
+ void init_members(void);
+ void copy_members(const ctlikelihood& app);
+ void free_members(void);
+ double evaluate(GModelPar& par, const double& value);
+
+ // Protected members
+ GObservations m_obs; //!< Observation container
+ GOptimizerLM m_opt; //!< Optimizer
+};
+
+
+/***********************************************************************//**
+ * @brief Return observation container
+ *
+ * @return Reference to observation container.
+ *
+ * Returns a reference to the observation container.
+ ***************************************************************************/
+inline
+const GObservations& ctlikelihood::obs(void) const
+{
+ return m_obs;
+}
+
+
+/***********************************************************************//**
+ * @brief Return optimizer
+ *
+ * @return Pointer to optimizer.
+ *
+ * Returns a pointer to the optimizer.
+ ***************************************************************************/
+inline
+const GOptimizer* ctlikelihood::opt(void) const
+{
+ return &m_opt;
+}
+
+#endif /* CTLIKELIHOOD_HPP */
diff --git a/src/support/ctool.cpp b/src/support/ctool.cpp
index f69018d..c37f159 100644
--- a/src/support/ctool.cpp
+++ b/src/support/ctool.cpp
@@ -1655,8 +1655,8 @@ GEnergies ctool::insert_energy_boundaries(const GEnergies& energies,
/***********************************************************************//**
* @brief Determine the counts cube layer usage.
*
- * @param[in] ebounds_cube Energy boundaries of the counts cube.
- * @param[in] ebounds_list Energy boundaries of the event list.
+ * @param[in] cube_ebounds Energy boundaries of the counts cube.
+ * @param[in] list_ebounds Energy boundaries of the event list.
* @return Vector of usage flags.
*
* Determines a vector of counts cube layer usage flags that signal whether
diff --git a/src/support/ctool.hpp b/src/support/ctool.hpp
index 4cadfca..f6ead66 100644
--- a/src/support/ctool.hpp
+++ b/src/support/ctool.hpp
@@ -37,7 +37,13 @@
/***********************************************************************//**
* @class ctool
*
- * @brief ctool base class
+ * @brief Base class for ctools
+ *
+ * This class is the generic base class for all ctools. It enforces
+ * implementation of the clear(), run() and save() methods by all ctools,
+ * and it implements the execute() method.
+ *
+ * The class also provides a number of helper methods for ctools.
***************************************************************************/
class ctool : public GApplication {
public:
diff --git a/test/test_ctulimit.py b/test/test_ctulimit.py
index 66705b9..9219c0f 100755
--- a/test/test_ctulimit.py
+++ b/test/test_ctulimit.py
@@ -103,9 +103,9 @@ class Test(test):
ulimit = ctools.ctulimit()
# Check that empty ctulimit tool holds zero upper limits
- self.test_value(ulimit.diff_ulimit(), 0.0, 1.0e-23,
+ self.test_value(ulimit.diff_ulimit(), 0.0, 1.0e-21,
'Check differential upper limit')
- self.test_value(ulimit.flux_ulimit(), 0.0, 1.0e-17,
+ self.test_value(ulimit.flux_ulimit(), 0.0, 1.0e-16,
'Check upper limit on photon flux')
self.test_value(ulimit.eflux_ulimit(), 0.0, 1.0e-16,
'Check upper limit on energy flux')
@@ -132,9 +132,9 @@ class Test(test):
ulimit.save()
# Check results
- self.test_value(ulimit.diff_ulimit(), 8.86803e-18, 1.0e-23,
+ self.test_value(ulimit.diff_ulimit(), 8.86803e-18, 1.0e-21,
'Check differential upper limit')
- self.test_value(ulimit.flux_ulimit(), 6.10509e-12, 1.0e-17,
+ self.test_value(ulimit.flux_ulimit(), 6.10509e-12, 1.0e-16,
'Check upper limit on photon flux')
self.test_value(ulimit.eflux_ulimit(), 2.75669e-11, 1.0e-16,
'Check upper limit on energy flux')
@@ -150,10 +150,10 @@ class Test(test):
cpy_ulimit = ulimit.copy()
# Check results of copy
- self.test_value(cpy_ulimit.diff_ulimit(), 8.86803e-18, 1.0e-23,
+ self.test_value(cpy_ulimit.diff_ulimit(), 8.86803e-18, 1.0e-21,
'Check differential upper limit')
- self.test_value(cpy_ulimit.flux_ulimit(), 6.10509e-12, 1.0e-17,
- 'ulimit upper limit on photon flux')
+ self.test_value(cpy_ulimit.flux_ulimit(), 6.10509e-12, 1.0e-16,
+ 'Check upper limit on photon flux')
self.test_value(cpy_ulimit.eflux_ulimit(), 2.75669e-11, 1.0e-16,
'Check upper limit on energy flux')
@@ -165,10 +165,10 @@ class Test(test):
cpy_ulimit.execute()
# Check results
- self.test_value(cpy_ulimit.diff_ulimit(), 8.86803e-18, 1.0e-23,
+ self.test_value(cpy_ulimit.diff_ulimit(), 8.86803e-18, 1.0e-21,
'Check differential upper limit')
- self.test_value(cpy_ulimit.flux_ulimit(), 6.10509e-12, 1.0e-17,
- 'ulimit upper limit on photon flux')
+ self.test_value(cpy_ulimit.flux_ulimit(), 6.10509e-12, 1.0e-16,
+ 'Check upper limit on photon flux')
self.test_value(cpy_ulimit.eflux_ulimit(), 2.75669e-11, 1.0e-16,
'Check upper limit on energy flux')
@@ -179,9 +179,9 @@ class Test(test):
cpy_ulimit.clear()
# Check that empty ctulimit tool holds zero upper limits
- self.test_value(cpy_ulimit.diff_ulimit(), 0.0, 1.0e-23,
+ self.test_value(cpy_ulimit.diff_ulimit(), 0.0, 1.0e-21,
'Check differential upper limit')
- self.test_value(cpy_ulimit.flux_ulimit(), 0.0, 1.0e-17,
+ self.test_value(cpy_ulimit.flux_ulimit(), 0.0, 1.0e-16,
'Check upper limit on photon flux')
self.test_value(cpy_ulimit.eflux_ulimit(), 0.0, 1.0e-16,
'Check upper limit on energy flux')
@@ -190,11 +190,7 @@ class Test(test):
# time. This should avoid the necessity to recompute the maximum
# likelihood
ulimit = ctools.ctulimit(obs)
- #ulimit['inobs'] = self._events
- #ulimit['inmodel'] = self._model
ulimit['srcname'] = 'Crab'
- #ulimit['caldb'] = self._caldb
- #ulimit['irf'] = self._irf
ulimit['logfile'] = 'ctulimit_py3.log'
ulimit['chatter'] = 4
@@ -203,10 +199,10 @@ class Test(test):
ulimit.execute()
# Check results
- self.test_value(ulimit.diff_ulimit(), 8.86803e-18, 1.0e-23,
+ self.test_value(ulimit.diff_ulimit(), 8.86803e-18, 1.0e-21,
'Check differential upper limit')
- self.test_value(ulimit.flux_ulimit(), 6.10509e-12, 1.0e-17,
- 'ulimit upper limit on photon flux')
+ self.test_value(ulimit.flux_ulimit(), 6.10509e-12, 1.0e-16,
+ 'Check upper limit on photon flux')
self.test_value(ulimit.eflux_ulimit(), 2.75669e-11, 1.0e-16,
'Check upper limit on energy flux')