34 #include "GOptimizer.hpp"
42 #define G_ERR_BISECTION "cterror::error_bisection(double&, double&)"
202 this->GApplication::clear();
236 log_header1(TERSE,
"Compute best-fit likelihood");
244 GOptimizerLM best_opt =
m_opt;
248 log_value(NORMAL,
"Maximum log likelihood", gammalib::str(
m_best_logL,3));
252 GModels models_best =
m_obs.models();
258 int npars = model->size();
261 for (
int i = 0; i < npars; ++i) {
264 if (model->at(i).is_fixed()) {
269 m_obs.models(models_best);
272 GModels& current_models =
const_cast<GModels&
>(
m_obs.models());
279 double parmin = std::max(
m_model_par->factor_min(),
281 double parmax = std::min(
m_model_par->factor_max(),
285 log_header1(TERSE,
"Compute error for source \""+
m_srcname+
"\""
287 log_value(NORMAL,
"Confidence level",
289 log_value(NORMAL,
"Log-likelihood difference",
m_dlogL);
290 log_value(NORMAL,
"Initial factor range",
291 "["+gammalib::str(parmin)+
", "+gammalib::str(parmax)+
"]");
298 double error = 0.5 * (value_hi - value_lo);
299 double error_neg =
m_value - value_lo;
300 double error_pos = value_hi -
m_value;
301 double error_value = std::abs(error*
m_model_par->scale());
302 double error_value_neg = std::abs(error_neg*
m_model_par->scale());
303 double error_value_pos = std::abs(error_pos*
m_model_par->scale());
307 log_value(NORMAL,
"Lower parameter factor", value_lo);
308 log_value(NORMAL,
"Upper parameter factor", value_hi);
309 log_value(NORMAL,
"Error from curvature",
311 log_value(NORMAL,
"Error from profile",
312 gammalib::str(error_value) + unit);
313 log_value(NORMAL,
"Negative profile error",
314 gammalib::str(error_value_neg) + unit);
315 log_value(NORMAL,
"Positive profile error",
316 gammalib::str(error_value_pos) + unit);
319 model->at(i).factor_error(error);
324 m_obs.models(models_best);
345 log_header1(TERSE,
"Save results");
348 if ((*
this)[
"outmodel"].is_valid()) {
354 log_value(NORMAL,
"Model definition file",
m_outmodel.url());
363 log_value(NORMAL,
"Model definition file",
"NONE");
463 double sigma = gammalib::erfinv(
m_confidence) * gammalib::sqrt_two;
467 m_opt.eps((*
this)[
"like_accuracy"].real());
468 m_opt.max_iter((*
this)[
"max_iter"].integer());
471 m_tol = (*this)[
"tol"].real();
473 m_chatter =
static_cast<GChatter
>((*this)[
"chatter"].integer());
477 (*this)[
"outmodel"].query();
482 int nthreads = (*this)[
"nthreads"].integer();
484 omp_set_num_threads(nthreads);
507 double wrk_min = min;
508 double wrk_max = max;
514 double mid = (wrk_min + wrk_max) / 2.0;
522 std::string msg =
"The \""+
m_model_par->name()+
"\" parameter "
523 "minimum has been reached during error "
524 "calculation. To obtain accurate errors, "
525 "consider setting the minimum parameter "
526 "value to a lower value, and re-run "
528 log_string(TERSE, msg);
532 std::string msg =
"The \""+
m_model_par->name()+
"\" parameter "
533 "maximum has been reached during error "
534 "calculation. To obtain accurate errors, "
535 "consider setting the maximum parameter "
536 "value to a higher value, and re-run "
538 log_string(TERSE, msg);
542 std::string msg =
"The maximum number of "+
544 "been reached. Please increase the "
545 "\"max_iter\" parameter, and re-run "
552 mid = (wrk_min + wrk_max) / 2.0;
558 log_value(EXPLICIT,
" Iteration "+gammalib::str(iter),
559 "["+gammalib::str(wrk_min)+
", "+gammalib::str(wrk_max)+
"]");
562 if (std::abs(eval_mid) <
m_tol) {
567 if (std::abs(wrk_max-wrk_min) < 1.0e-6) {
575 if (eval_mid > 0.0) {
578 else if (eval_mid < 0.0) {
587 if (eval_mid > 0.0) {
590 else if (eval_mid < 0.0) {
void save(void)
Save model.
Parameter error calculation tool.
void free_members(void)
Delete class members.
double m_tol
Tolerance for limit determination.
ctlikelihood & operator=(const ctlikelihood &app)
Assignment operator.
void get_parameters(void)
Get application parameters.
double evaluate(GModelPar &par, const double &value)
Evaluates the log-likelihood function.
GChatter m_chatter
Chattiness.
double m_dlogL
Likelihood difference for upper limit computation.
double error_bisection(const double &min, const double &max)
Calculate error using a bisection method.
void init_members(void)
Initialise class members.
virtual ~cterror(void)
Destructor.
void set_obs_statistic(const std::string &statistic)
Set fit statistic for CTA observations.
double m_value
Parameter value.
std::string m_srcname
Name of source.
void process(void)
Compute parameter errors using a likelihood profile method.
cterror(void)
Void constructor.
void clear(void)
Clear cterror tool.
Parameter error calculation tool interface definition.
GModelPar * m_model_par
Pointer to model parameter.
GFilename m_outmodel
Output model XML file.
void free_members(void)
Delete class members.
void free_members(void)
Delete class members.
bool m_apply_edisp
Apply energy dispersion?
void init_members(void)
Initialise class members.
void copy_members(const cterror &app)
Copy class members.
void init_members(void)
Initialise class members.
cterror & operator=(const cterror &app)
Assignment operator.
GObservations m_obs
Observation container.
double m_confidence
Confidence level.
int m_max_iter
Maximum number of iterations.
double m_best_logL
Best fit log likelihood of given model.
Base class for likelihood tools.
GOptimizerLM m_opt
Optimizer.