Add accept_dec parameter to likelihood fitting tools and use best observation in SRCLIX methods

Jürgen Knödlseder
1 parent 3d26c9ce
Showing 10 changed files with 116 additions and 35 deletions
doc/source/users/reference_manual/comlixfit.rst
doc/source/users/reference_manual/comlixmap.rst
doc/source/users/reference_manual/ctlike.rst
modules/comscripts/comlixfit.par
modules/comscripts/comlixfit.py
modules/comscripts/comlixmap.par
modules/comscripts/comlixmap.py
src/ctlike/ctlike.cpp
src/ctlike/ctlike.hpp
src/ctlike/ctlike.par
@@ -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?
@@ -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?
@@ -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.
@@ -28,6 +28,7 @@ outmodel,  f, a, results.xml,,, &quot;Output model definition file&quot;
 #==================
 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"
@@ -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)
@@ -37,6 +37,7 @@ binsz,    r, a, 1.0,0.01,10.0, &quot;Image scale (in degrees/pixel)&quot;
 #==================
 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"
@@ -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()):
 /***************************************************************************
  *                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&amp; 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();
@@ -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?
@@ -35,6 +35,7 @@ statistic,       s, h, DEFAULT,DEFAULT|CSTAT|WSTAT|CHI2,, &quot;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?"