From 63cb462cbedf3fb028274518082014fc59f5c1e0 Mon Sep 17 00:00:00 2001 From: Jürgen Knödlseder Date: Mon, 14 Feb 2022 12:00:37 +0100 Subject: [PATCH] Add accept_dec parameter to likelihood fitting tools and use best observation in SRCLIX methods --- doc/source/users/reference_manual/comlixfit.rst | 9 ++++++++- doc/source/users/reference_manual/comlixmap.rst | 9 ++++++++- doc/source/users/reference_manual/ctlike.rst | 6 ++++++ modules/comscripts/comlixfit.par | 1 + modules/comscripts/comlixfit.py | 34 +++++++++++++++++++++++++++------- modules/comscripts/comlixmap.par | 1 + modules/comscripts/comlixmap.py | 79 ++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++------------------- src/ctlike/ctlike.cpp | 10 ++++------ src/ctlike/ctlike.hpp | 1 - src/ctlike/ctlike.par | 1 + 10 files changed, 116 insertions(+), 35 deletions(-) diff --git a/doc/source/users/reference_manual/comlixfit.rst b/doc/source/users/reference_manual/comlixfit.rst index 23df08f..a760323 100755 --- a/doc/source/users/reference_manual/comlixfit.rst +++ b/doc/source/users/reference_manual/comlixfit.rst @@ -38,7 +38,14 @@ General parameters Maximum number of SRCLIX iterations. ``(like_accuracy = 0.05) [real]`` - Absolute accuracy of maximum likelihood value. + Absolute accuracy of maximum likelihood value. The SRCLIX iterations terminate + if the log-likelihood improvement is below this accuracy. + +``(accept_dec = 0.0) [real]`` + Maximum accepted log-likelihood decrease. Setting this parameter to a positive + value with allow some decrease of the log-likelihood. This may help to get the + algorithm out of a local minimum. Use this parameter with case since it may + lead to a solution that is not the maximum likelihood solution. ``(fix_spat_for_ts = no) [boolean]`` Fix spatial parameters for TS computation? diff --git a/doc/source/users/reference_manual/comlixmap.rst b/doc/source/users/reference_manual/comlixmap.rst index 63a14d6..01b42e2 100755 --- a/doc/source/users/reference_manual/comlixmap.rst +++ b/doc/source/users/reference_manual/comlixmap.rst @@ -31,7 +31,14 @@ General parameters Maximum number of SRCLIX iterations. ``(like_accuracy = 0.05) [real]`` - Absolute accuracy of maximum likelihood value. + Absolute accuracy of maximum likelihood value. The SRCLIX iterations terminate + if the log-likelihood improvement is below this accuracy. + +``(accept_dec = 0.0) [real]`` + Maximum accepted log-likelihood decrease. Setting this parameter to a positive + value with allow some decrease of the log-likelihood. This may help to get the + algorithm out of a local minimum. Use this parameter with case since it may + lead to a solution that is not the maximum likelihood solution. ``(fix_spat_for_ts = no) [boolean]`` Fix spatial parameters for TS computation? diff --git a/doc/source/users/reference_manual/ctlike.rst b/doc/source/users/reference_manual/ctlike.rst index 4eecf3f..fb0d3ec 100755 --- a/doc/source/users/reference_manual/ctlike.rst +++ b/doc/source/users/reference_manual/ctlike.rst @@ -128,6 +128,12 @@ General parameters log likelihood value. Converserly, decreasing the value will result in less iterations at the expense of a less accurate maximum likelihood value. +``(accept_dec = 0.0) [real]`` + Maximum accepted log-likelihood decrease. Setting this parameter to a positive + value with allow some decrease of the log-likelihood. This may help to get the + algorithm out of a local minimum. Use this parameter with case since it may + lead to a solution that is not the maximum likelihood solution. + ``(max_iter = 50) [integer]`` Maximum number of fit iterations. diff --git a/modules/comscripts/comlixfit.par b/modules/comscripts/comlixfit.par index aa193ad..ef773a3 100644 --- a/modules/comscripts/comlixfit.par +++ b/modules/comscripts/comlixfit.par @@ -28,6 +28,7 @@ outmodel, f, a, results.xml,,, "Output model definition file" #================== max_iter, i, h, 50,1,1000, "Maximum number of SRCLIX iterations" like_accuracy, r, h, 0.05,1.0e-30,1.0, "Absolute accuracy of maximum likelihood value" +accept_dec, r, h, 0.0,0.0,10.0, "Maximum accepted log-likelihood decrease" fix_spat_for_ts, b, h, no,,, "Fix spatial parameters for TS computation?" fix_ext_for_initial, b, h, no,,, "Fix extension parameters for initial fit iteration?" bkgmethod, s, a, BGDLIXE,PHINOR|BGDLIXA|BGDLIXE,, "Method for background computation" diff --git a/modules/comscripts/comlixfit.py b/modules/comscripts/comlixfit.py index a7c6cf8..5fab2cb 100755 --- a/modules/comscripts/comlixfit.py +++ b/modules/comscripts/comlixfit.py @@ -65,6 +65,7 @@ class comlixfit(ctools.cslikelihood): # Query parameters self['max_iter'].integer() self['like_accuracy'].real() + self['accept_dec'].real() self['fix_spat_for_ts'].boolean() self['fix_ext_for_initial'].boolean() @@ -265,6 +266,11 @@ class comlixfit(ctools.cslikelihood): self._log_header1(gammalib.NORMAL, 'Iterative maximum likelihood model fitting') + # Initialise optimiser + self.opt().max_iter(100) + self.opt().eps(0.005) + self.opt().accept_dec(self['accept_dec'].real()) + # Loop over iterations for iter in range(niter): @@ -274,6 +280,10 @@ class comlixfit(ctools.cslikelihood): # Update observations self._update_obs() + # If first iteration then initialise best observation + if iter == 0: + best_obs = self.obs().copy() + # Optionally fix or free fixed spatial extension parameters if self['fix_ext_for_initial'].boolean(): if iter == 0: @@ -284,21 +294,26 @@ class comlixfit(ctools.cslikelihood): # Fit model self.obs().optimize(self.opt()) - # Compute logL difference after first iteration + # Compute logL difference if iter > 0: delta = logL - self.opt().value() + else: + logL = self.opt().value() - # Store maximum likelihood value - logL = self.opt().value() + # If log-likelihood improved then update best observations and + # last log-likelihood value + if delta >= 0.0: + best_obs = self.obs().copy() + logL = self.opt().value() # Increment iterations counter self._iter += self.opt().iter() # Log maximum likelihood if iter == 0: - result = '%.5f' % (logL) + result = '%.5f' % (self.opt().value()) else: - result = '%.5f (%.5f)' % (logL, delta) + result = '%.5f (%.5f)' % (self.opt().value(), delta) self._log_value(gammalib.NORMAL, 'logL after iteration %d' % (iter+1), result) @@ -314,12 +329,17 @@ class comlixfit(ctools.cslikelihood): self._log_header1(gammalib.EXPLICIT, 'Final maximum likelihood model fitting') + # Use best observations for maximum likelihood fitting. This avoids + # that observations that have a worse log-likelihood than the best + # value are used for final model fitting. + self.obs(best_obs) + # Do final model fit self._final_model_fit() # Compute logL difference and store maximum likelihood value - delta = logL - self.opt().value() - logL = self.opt().value() + delta = logL - self.opt().value() + logL = self.opt().value() # Log final maximum likelihood result = '%.5f (%.5f)' % (logL, delta) diff --git a/modules/comscripts/comlixmap.par b/modules/comscripts/comlixmap.par index 167e09d..4c4499f 100644 --- a/modules/comscripts/comlixmap.par +++ b/modules/comscripts/comlixmap.par @@ -37,6 +37,7 @@ binsz, r, a, 1.0,0.01,10.0, "Image scale (in degrees/pixel)" #================== max_iter, i, h, 50,1,1000, "Maximum number of SRCLIX iterations" like_accuracy, r, h, 0.05,1.0e-30,1.0, "Absolute accuracy of maximum likelihood value" +accept_dec, r, h, 0.0,0.0,10.0, "Maximum accepted log-likelihood decrease" fix_spat, b, h, yes,,, "Fix spatial parameters for TS map computation?" bkgmethod, s, a, BGDLIXE,PHINOR|BGDLIXA|BGDLIXE,, "Method for background computation" nrunav, i, h, 3,,, "Number of Chi/Psi bins used for running average" diff --git a/modules/comscripts/comlixmap.py b/modules/comscripts/comlixmap.py index e423d5d..1779ac7 100755 --- a/modules/comscripts/comlixmap.py +++ b/modules/comscripts/comlixmap.py @@ -69,6 +69,7 @@ class comlixmap(ctools.cslikelihood): # Query parameters self['like_accuracy'].real() self['max_iter'].integer() + self['accept_dec'].real() self['fix_spat'].boolean() self['coordsys'].string() self['proj'].string() @@ -210,6 +211,9 @@ class comlixmap(ctools.cslikelihood): ts : float TS of test source """ + # Get list of flux normalisation parameters + pars = ['Prefactor', 'Normalization', 'EnergyFlux', 'PhotonFlux'] + # Loop over test source names for srcname in self._srcnames: @@ -245,21 +249,31 @@ class comlixmap(ctools.cslikelihood): # Update observations self._update_obs() + # If first iteration then initialise best observation + if iter == 0: + best_obs = self.obs().copy() + # Fit model self.obs().optimize(self.opt()) - # Compute logL difference after first iteration + # Compute log-likelihood difference or initialise + # log-likelihood value if iter > 0: delta = logL - self.opt().value() + else: + logL = self.opt().value() - # Store maximum likelihood value - logL = self.opt().value() + # If log-likelihood improved then update best observations and + # last log-likelihood value + if delta >= 0.0: + best_obs = self.obs().copy() + logL = self.opt().value() # Log maximum likelihood if iter == 0: - result = '%.5f' % (logL) + result = '%.5f' % (self.opt().value()) else: - result = '%.5f (%.5f)' % (logL, delta) + result = '%.5f (%.5f)' % (self.opt().value(), delta) self._log_value(gammalib.NORMAL, 'logL after iteration %d' % (iter+1), result) @@ -268,13 +282,37 @@ class comlixmap(ctools.cslikelihood): if delta < eps: break - # Compute TS in final model fit - ts = self._final_model_fit() + # Use best observations for maximum likelihood fitting. This avoids + # that observations that have a worse log-likelihood than the best + # value are used for final model fitting. + self.obs(best_obs) + + # Compute log-likelihood and TS for best observations + value = self._final_model_fit() - # Log TS - for i, srcname in enumerate(self._srcnames): - key = 'TS %s' % srcname - self._log_value(gammalib.NORMAL, key, ts[i]) + # Compute log-likelihood difference + delta = logL - value + + # Log final maximum likelihood + result = '%.5f (%.5f)' % (value, delta) + self._log_value(gammalib.NORMAL, 'logL after final iteration', + result) + + # Collect and log TS and source normalisation + ts = [] + for srcname in self._srcnames: + source = self.obs().models()[srcname] + ts_value = float(source.ts()) + key = 'TS %s' % srcname + value = '%.3f' % ts_value + ts.append(ts_value) + for par in pars: + if source.has_par(par): + value += ' (%e +/- %e %s)' % (source[par].value(), + source[par].error(), + source[par].unit()) + break + self._log_value(gammalib.NORMAL, key, value) # Return TS return ts @@ -285,8 +323,8 @@ class comlixmap(ctools.cslikelihood): Returns ------- - ts : float - TS of test source + logL : float + Log-likelihood of final model fit """ # Create instance of model fitting tool like = ctools.ctlike(self.obs()) @@ -295,13 +333,11 @@ class comlixmap(ctools.cslikelihood): # Run ctlike like.run() - # Recover TS for all test sources - ts = [] - for srcname in self._srcnames: - ts.append(like.obs().models()[srcname].ts()) + # Recover results + self.obs(like.obs()) - # Return TS - return ts + # Return log-likelihood value + return (like.opt().value()) # Public methods @@ -357,6 +393,11 @@ class comlixmap(ctools.cslikelihood): # Get a copy of the initial models models = self.obs().models().copy() + # Initialise optimiser + self.opt().max_iter(100) + self.opt().eps(0.005) + self.opt().accept_dec(self['accept_dec'].real()) + # Loop over grid positions for ipix in range(self._maps[0].npix()): diff --git a/src/ctlike/ctlike.cpp b/src/ctlike/ctlike.cpp index 190b010..11243fc 100755 --- a/src/ctlike/ctlike.cpp +++ b/src/ctlike/ctlike.cpp @@ -1,7 +1,7 @@ /*************************************************************************** * ctlike - Maximum likelihood fitting tool * * ----------------------------------------------------------------------- * - * copyright (C) 2010-2020 by Juergen Knoedlseder * + * copyright (C) 2010-2022 by Juergen Knoedlseder * * ----------------------------------------------------------------------- * * * * This program is free software: you can redistribute it and/or modify * @@ -425,7 +425,6 @@ void ctlike::init_members(void) m_outmodel.clear(); m_outcovmat.clear(); m_max_iter = 50; - m_like_accuracy = 0.005; m_refit = false; m_refit_if_failed = true; m_apply_edisp = false; @@ -455,7 +454,6 @@ void ctlike::copy_members(const ctlike& app) m_outmodel = app.m_outmodel; m_outcovmat = app.m_outcovmat; m_max_iter = app.m_max_iter; - m_like_accuracy = app.m_like_accuracy; m_refit = app.m_refit; m_refit_if_failed = app.m_refit_if_failed; m_apply_edisp = app.m_apply_edisp; @@ -507,10 +505,10 @@ void ctlike::get_parameters(void) } // endif: no models were associated with observations // Set optimizer characteristics from user parameters - m_max_iter = (*this)["max_iter"].integer(); - m_like_accuracy = (*this)["like_accuracy"].real(); - m_opt.eps(m_like_accuracy); + m_max_iter = (*this)["max_iter"].integer(); m_opt.max_iter(m_max_iter); + m_opt.eps((*this)["like_accuracy"].real()); + m_opt.accept_dec((*this)["accept_dec"].real()); // Get other parameters m_refit = (*this)["refit"].boolean(); diff --git a/src/ctlike/ctlike.hpp b/src/ctlike/ctlike.hpp index bbd3ae1..180f4c4 100755 --- a/src/ctlike/ctlike.hpp +++ b/src/ctlike/ctlike.hpp @@ -76,7 +76,6 @@ protected: GFilename m_outmodel; //!< Source model output XML file name GFilename m_outcovmat; //!< Covariance matrix output file name int m_max_iter; //!< Maximum number of iterations - double m_like_accuracy; //!< Log-likelihood accuracy bool m_refit; //!< Refitting? bool m_refit_if_failed; //!< Refitting in case of failure? bool m_apply_edisp; //!< Apply energy dispersion? diff --git a/src/ctlike/ctlike.par b/src/ctlike/ctlike.par index 854afca..3a9097d 100755 --- a/src/ctlike/ctlike.par +++ b/src/ctlike/ctlike.par @@ -35,6 +35,7 @@ statistic, s, h, DEFAULT,DEFAULT|CSTAT|WSTAT|CHI2,, "Optimization statisti refit, b, h, no,,, "Do refitting?" refit_if_failed, b, h, yes,,, "Refit if fit failed?" like_accuracy, r, h, 0.005,1.0e-30,1.0, "Absolute accuracy of maximum likelihood value" +accept_dec, r, h, 0.0,0.0,10.0, "Maximum accepted log-likelihood decrease" max_iter, i, h, 50,1,1000, "Maximum number of fit iterations" fix_spat_for_ts, b, h, no,,, "Fix spatial parameters for TS computation?" -- libgit2 0.21.4