Karim Sabri / ctools

Commit d223acf4b2a1b2b5323422da19634445462ac2ba

Authored by Jürgen Knödlseder
1 parent 4109dc8f

Re-compute best log-likelihood value for fixed value and exit iterations if para… 

…meter interval becomes small
Inline Side-by-side
Showing 1 changed file with 31 additions and 3 deletions
src/ctulimit/ctulimit.cpp
  View file @d223acf
1 /*************************************************************************** 1 /***************************************************************************
2 * ctulimit - Upper limit calculation tool * 2 * ctulimit - Upper limit calculation tool *
3 * ----------------------------------------------------------------------- * 3 * ----------------------------------------------------------------------- *
4 - * copyright (C) 2015-2020 by Michael Mayer * 4 + * copyright (C) 2015-2021 by Michael Mayer *
5 * ----------------------------------------------------------------------- * 5 * ----------------------------------------------------------------------- *
6 * * 6 * *
7 * This program is free software: you can redistribute it and/or modify * 7 * This program is free software: you can redistribute it and/or modify *
@@ -268,7 +268,8 @@ void ctulimit::run(void) @@ -268,7 +268,8 @@ void ctulimit::run(void)
268 268
269 } // endif: likelihood was zero 269 } // endif: likelihood was zero
270 270
271 - // Get parameter brackets 271 + // Get parameter brackets. Note that this method recomputes m_best_logL
  272 + // to make sure that the value is coherent for the upper limit search.
272 double parmin; 273 double parmin;
273 double parmax; 274 double parmax;
274 get_parameter_brackets(parmin, parmax); 275 get_parameter_brackets(parmin, parmax);
@@ -280,6 +281,7 @@ void ctulimit::run(void) @@ -280,6 +281,7 @@ void ctulimit::run(void)
280 log_value(NORMAL, "Model name", m_skymodel->name()); 281 log_value(NORMAL, "Model name", m_skymodel->name());
281 log_value(NORMAL, "Parameter name", m_model_par->name()); 282 log_value(NORMAL, "Parameter name", m_model_par->name());
282 log_value(NORMAL, "Confidence level", gammalib::str(m_confidence*100.0)+"%"); 283 log_value(NORMAL, "Confidence level", gammalib::str(m_confidence*100.0)+"%");
  284 + log_value(NORMAL, "Best log-likelihood value", m_best_logL);
283 log_value(NORMAL, "Log-likelihood difference", m_dlogL); 285 log_value(NORMAL, "Log-likelihood difference", m_dlogL);
284 log_value(NORMAL, "Initial parameter range", 286 log_value(NORMAL, "Initial parameter range",
285 "["+gammalib::str(parmin)+", "+gammalib::str(parmax)+"]"); 287 "["+gammalib::str(parmin)+", "+gammalib::str(parmax)+"]");
@@ -604,6 +606,13 @@ void ctulimit::get_parameter_brackets(double& parmin, double& parmax) @@ -604,6 +606,13 @@ void ctulimit::get_parameter_brackets(double& parmin, double& parmax)
604 double value = m_model_par->factor_value(); 606 double value = m_model_par->factor_value();
605 double error = m_model_par->factor_error(); 607 double error = m_model_par->factor_error();
606 608
  609 + // Re-compute best log-likelihood for value fixed at optimum
  610 + m_best_logL = evaluate(*m_model_par, value);
  611 +
  612 + // Write best log-likelihood
  613 + log_value(NORMAL, "Best log-likelihood",
  614 + "logL("+gammalib::str(value)+")="+gammalib::str(m_best_logL));
  615 +
607 // If parameter error is zero then take parameter value as error 616 // If parameter error is zero then take parameter value as error
608 if (error == 0.0) { 617 if (error == 0.0) {
609 error = value; 618 error = value;
@@ -671,7 +680,12 @@ void ctulimit::get_parameter_brackets(double& parmin, double& parmax) @@ -671,7 +680,12 @@ void ctulimit::get_parameter_brackets(double& parmin, double& parmax)
671 * Exceeded maximum number of iterations in the upper limit 680 * Exceeded maximum number of iterations in the upper limit
672 * search. 681 * search.
673 * 682 *
674 - * Calculates the upper limit using a bisection method. 683 + * Calculates the upper limit using a bisection method. The bisection is
  684 + * stopped if the log-likelihood value is within the tolerance of the target
  685 + * value. The bisection is also stopped putting a warning in the log file if
  686 + * the parameter interval becomes smaller than a tolerance of 1.0e-6. The
  687 + * method stops with an exception if the maximum number of iterations is
  688 + * exhausted.
675 ***************************************************************************/ 689 ***************************************************************************/
676 void ctulimit::ulimit_bisection(const double& parmin, const double& parmax) 690 void ctulimit::ulimit_bisection(const double& parmin, const double& parmax)
677 { 691 {
@@ -710,6 +724,20 @@ void ctulimit::ulimit_bisection(const double& parmin, const double& parmax) @@ -710,6 +724,20 @@ void ctulimit::ulimit_bisection(const double& parmin, const double& parmax)
710 break; 724 break;
711 } 725 }
712 726
  727 + // Check for vanishing parameter range
  728 + if (std::abs(wrk_min-wrk_max) < 1.0e-6) {
  729 + std::string msg = "*** WARNING: Parameter range "
  730 + "["+gammalib::str(wrk_min)+", "+
  731 + gammalib::str(wrk_max)+"] has reduced to a "
  732 + "narrow interval without reaching convergence!\n"
  733 + " The best log-likelihood value is "
  734 + "probably not an absolute but only a local "
  735 + "minimum.\n"
  736 + " Stop iterations.";
  737 + log_string(TERSE, msg);
  738 + break;
  739 + }
  740 +
713 // Change boundaries for further iteration 741 // Change boundaries for further iteration
714 if (eval_mid > 0.0) { 742 if (eval_mid > 0.0) {
715 wrk_max = mid; // logL too large, new range = [wrk_min, mid] 743 wrk_max = mid; // logL too large, new range = [wrk_min, mid]