Integrate analytical gradients for energy dispersion computation (#3333)

Jürgen Knödlseder
1 parent 81862a72
Showing 10 changed files with 235 additions and 78 deletions
ChangeLog
NEWS
README.md
doc/source/admin/release_history/1.8.rst
include/GIntegrals.hpp
include/GResponse.hpp
include/Makefile.am
src/numerics/GIntegrals.cpp
src/obs/GObservation.cpp
src/obs/GResponse.cpp
-2020-09-04
+2020-09-11
  
         * Version 1.8.0 released
         ========================
  
+        Integrate analytical gradients for energy dispersion computation (#3333)
         Add GFunctions abstract base class (#3328)
         Add GIntegrals class (#3328)
  
 New Features and Important Changes in GammaLib 1.8.0
  
-4 September 2020
+11 September 2020
  
  
 1. Introduction
@@ -76,7 +76,14 @@ of functions (#3328).
  
 10. Observation module
 ----------------------
-None
+The GResponse::edisp_kern integration kernel class that is used for the integration
+over true energies in case that energy dispersion should be considered now derives
+from GFunctions. The actual function value is returned in the first slot of the
+GNdarray returned by the GResponse::edisp_kern::eval() method, while the other slots
+contain the analytical parameter gradients. The GIntegrals class is now used in the
+GResponse::convolve() method to perform the true energy integration, integrating not
+only the function but also any available parameter gradients (#3333).
+
  
  
 11. Optimization module
 GammaLib information
 ====================
-* Version: 1.8.0.dev (4 September 2020)
+* Version: 1.8.0.dev (11 September 2020)
  
 [![Build Status](https://cta-jenkins.irap.omp.eu/buildStatus/icon?job=gammalib-integrate-os)](https://cta-jenkins.irap.omp.eu/job/gammalib-integrate-os/)
  
 .. _1.8:
  
-GammaLib 1.8 release (3 September 2020)
-=======================================
+GammaLib 1.8 release (11 September 2020)
+========================================
  
 Overview
 --------
@@ -21,6 +21,8 @@ Bug fixes
 Improvements
 ------------
  
+* [`3333 <https://cta-redmine.irap.omp.eu/issues/3333>`_] -
+  Integrate analytical gradients for energy dispersion computation
  
 New features
 ------------
@@ -96,7 +96,7 @@ protected:
                   const int&      n,
                   const int&      index,
                   const double&   x,
-                  GNdarray*       dy);
+                  double*         dy);
  
     // Protected data area
     GFunctions* m_kernels;   //!< Pointer to function kernels
@@ -33,6 +33,8 @@
 #include "GFunction.hpp"
 #include "GTime.hpp"
 #include "GResponseCache.hpp"
+#include "GFunctions.hpp"
+#include "GNdarray.hpp"
  
 /* __ Forward declarations _______________________________________________ */
 class GMatrix;
@@ -43,6 +45,7 @@ class GSkyDir;
 class GEnergy;
 class GEbounds;
 class GObservation;
+class GModelPar;
 class GModelSky;
 class GModelSpatialRadial;
 class GModelSpatialElliptical;
@@ -136,28 +139,25 @@ protected:
                                  const GObservation& obs) const;
  
     // Protected classes
-    class edisp_kern : public GFunction {
+    class edisp_kern : public GFunctions {
     public:
         edisp_kern(const GResponse*    parent,
                    const GObservation* obs,
                    const GModelSky*    model,
                    const GEvent*       event,
                    const GTime&        srcTime,
-                   const bool&         grad) :
-                   m_parent(parent),
-                   m_obs(obs),
-                   m_model(model),
-                   m_event(event),
-                   m_srcTime(srcTime),
-                   m_grad(grad) { }
-        double eval(const double& etrue);
+                   const bool&         grad);
+        const GNdarray& array(void) const;
+        GNdarray        eval(const double& etrue);
     protected:
-        const GResponse*    m_parent;  //!< Response
-        const GObservation* m_obs;     //!< Observation
-        const GModelSky*    m_model;   //!< Sky model
-        const GEvent*       m_event;   //!< Event
-        GTime               m_srcTime; //!< True arrival time
-        bool                m_grad;    //!< Gradient flag
+        const GResponse*        m_parent;  //!< Response
+        const GObservation*     m_obs;     //!< Observation
+        const GModelSky*        m_model;   //!< Sky model
+        const GEvent*           m_event;   //!< Event
+        GNdarray                m_array;   //!< Array of values and gradients
+        std::vector<GModelPar*> m_pars;    //!< Parameter pointers
+        GTime                   m_srcTime; //!< True arrival time
+        bool                    m_grad;    //!< Gradient flag
     };
     class irf_radial_kern_theta : public GFunction {
     public:
@@ -41,8 +41,10 @@ pkginclude_HEADERS = GammaLib.hpp \
                      GMatrixSparse.hpp \
                      GMatrixSymmetric.hpp \
                      GIntegral.hpp \
+                     GIntegrals.hpp \
                      GDerivative.hpp \
                      GFunction.hpp \
+                     GFunctions.hpp \
                      GMath.hpp \
                      GNdarray.hpp \
                      GFft.hpp \
@@ -303,11 +303,12 @@ GNdarray GIntegrals::romberg(const double&amp; a, const double&amp; b, const int&amp; order)
         }
  
         // Allocate temporal storage
-        GNdarray* s = new GNdarray[max_iter+2];
         double*   h = new double[max_iter+2];
+        GNdarray* s = new GNdarray[max_iter+2];
  
-        // Initialise step size
+        // Initialise step size and initial result
         h[1] = 1.0;
+        s[0] = GNdarray(m_kernels->array().shape());
  
         // Iterative loop
         for (m_iter = 1; m_iter <= max_iter; ++m_iter) {
@@ -333,9 +334,10 @@ GNdarray GIntegrals::romberg(const double&amp; a, const double&amp; b, const int&amp; order)
  
                 // Compute result using polynom interpolation
                 for (int i = 0; i < result.size(); ++i) {
-                    result(i) = polint(&h[m_iter-order],
-                                       &s[m_iter-order],
-                                       order, i, 0.0, &dss);
+                    double dy = 0;
+                    result(i) = polint(&h[m_iter-order], &s[m_iter-order],
+                                       order, i, 0.0, &dy);
+                    dss(i)    = dy;
                 }
  
                 // If a fixed number of iterations has been requested and if
@@ -369,8 +371,8 @@ GNdarray GIntegrals::romberg(const double&amp; a, const double&amp; b, const int&amp; order)
                     m_abserr     = abs(dss);
                     m_relerr     = m_abserr;
                     for (int i = 0; i < result.size(); ++i) {
-                        double abs_result = abs(result(i));
-                        if (abs_result != 0.0) {
+                        double abs_result = std::abs(result(i));
+                        if (abs_result > 0.0) {
                             m_relerr(i) /= abs_result;
                             m_has_relerr = true;
                         }
@@ -388,6 +390,10 @@ GNdarray GIntegrals::romberg(const double&amp; a, const double&amp; b, const int&amp; order)
  
         } // endfor: iterative loop
  
+        // Delete temporal storage
+        delete [] s;
+        delete [] h;
+
         // Set status and optionally dump warning
         if (!converged) {
             m_isvalid = false;
@@ -405,10 +411,6 @@ GNdarray GIntegrals::romberg(const double&amp; a, const double&amp; b, const int&amp; order)
             }
         }
  
-        // Delete temporal storage
-        delete [] s;
-        delete [] h;
-
     } // endif: integration range was valid
  
     // Return result
@@ -478,7 +480,7 @@ GNdarray GIntegrals::trapzd(const double&amp; a,
             m_calls += 2;
  
             // Compute result
-            result = 0.5*(b-a)*(y_a + y_b);
+            result = 0.5 * (b-a) * (y_a + y_b);
  
         } // endif: only a single step was requested
  
@@ -526,22 +528,20 @@ GNdarray GIntegrals::trapzd(const double&amp; a,
             }
  
             // Sum up values
-            double x   = a + 0.5*del;
+            double   x = a + 0.5 * del;
             GNdarray sum(result.shape());
-            for (int j = 0; j < it; ++j, x+=del) {
+            for (int j = 0; j < it; ++j, x += del) {
  
-                // Evaluate integrand
-                GNdarray y = m_kernels->eval(x);
+                // Evaluate and add integrand
+                sum += m_kernels->eval(x);
                 m_calls++;
  
-                // Add integrand
-                sum += y;
-
             } // endfor: looped over steps
  
             // Set result
-            result = 0.5*(result + (b-a)*sum/tnm);
-        }
+            result = 0.5 * (result + del * sum);
+
+        } // endelse: Case B
  
     } // endelse: trapeziodal rule was applied
  
@@ -708,7 +708,7 @@ double GIntegrals::polint(const double*   xa,
                           const int&      n,
                           const int&      index,
                           const double&   x,
-                          GNdarray*       dy)
+                          double*         dy)
 {
     // Initialise result
     double y = 0.0;
@@ -728,8 +728,9 @@ double GIntegrals::polint(const double*   xa,
             ns  = i;
             dif = dift;
         }
-        c[i] = ya[i+1](index);
-        d[i] = ya[i+1](index);
+        double value = ya[i+1](index);
+        c[i]         = value;
+        d[i]         = value;
     }
  
     // Get initial approximation to y
@@ -760,10 +761,10 @@ double GIntegrals::polint(const double*   xa,
         }
  
         // Compute y correction
-        (*dy)(index) = (2*(ns+1) < (n-m)) ? c[ns+1] : d[ns--];
+        *dy = (2*(ns+1) < (n-m)) ? c[ns+1] : d[ns--];
  
         // Update y
-        y += (*dy)(index);
+        y += *dy;
  
     } // endfor: looped over columns of tableau
  
@@ -271,9 +271,12 @@ double GObservation::model(const GModels&amp; models,
     int    igrad     = 0;      // Reset gradient counter
     int    grad_size = 0;      // Reset gradient size
  
+    // Set gradient usage flag
+    bool use_grad = (gradient != NULL);
+
     // If gradient is available then reset gradient vector elements to 0
     // and determine vector size
-    if (gradient != NULL) {
+    if (use_grad) {
         (*gradient) = 0.0;
         grad_size   = gradient->size();
     }
@@ -289,19 +292,11 @@ double GObservation::model(const GModels&amp; models,
             // observation identifier
             if (mptr->is_valid(instrument(), id())) {
  
-                // Compute value and add to model. If energy dispersion
-                // is used, don't compute model gradients as we cannot
-                // use them. This is somehow a kluge, but makes the
-                // code faster
-                model += mptr->eval(event, *this, !response()->use_edisp());
-
-                // Optionally determine model gradients. If the model has a
-                // gradient then use it, unless we have energy dispersion.
-                // For energy dispersion, no gradients are available as we
-                // have not implemented code that integrates the gradients
-                // over the energy dispersion. Here it's simpler to just use
-                // numerical gradients.
-                if (gradient != NULL) {
+                // Compute value and add to model
+                model += mptr->eval(event, *this, use_grad);
+
+                // Optionally determine model gradients
+                if (use_grad) {
  
                     // Make sure that we have a slot for the gradient
                     #if defined(G_RANGE_CHECK)
@@ -334,7 +329,7 @@ double GObservation::model(const GModels&amp; models,
  
                             // Set gradient
                             if (par.is_free()) {
-                                if (par.has_grad() && !response()->use_edisp()) {
+                                if (par.has_grad()) {
                                     (*gradient)[igrad+ipar] = par.factor_gradient();
                                 }
                                 else {
@@ -360,7 +355,7 @@ double GObservation::model(const GModels&amp; models,
  
                             // Set gradient
                             if (par.is_free()) {
-                                if (par.has_grad() && !response()->use_edisp()) {
+                                if (par.has_grad()) {
                                     (*gradient)[igrad+ipar] = par.factor_gradient();
                                 }
                                 else {
@@ -512,7 +507,7 @@ double GObservation::npred(const GModels&amp; models, GVector* gradient) const
  *                  S(p,E,t) \times R(p',E',t'|p,E,t) \, dp \, dE \, dt
  * \f]
  *
- * where                         
+ * where
  * \f$S(p,E,t)\f$ is the source model,
  * \f$R(p',E',t'|p,E,t)\f$ is the instrument response function,
  * \f$p'\f$ is the measured photon direction,
@@ -36,6 +36,7 @@
 #include "GVector.hpp"
 #include "GMatrix.hpp"
 #include "GIntegral.hpp"
+#include "GIntegrals.hpp"
 #include "GResponse.hpp"
 #include "GEvent.hpp"
 #include "GPhoton.hpp"
@@ -165,7 +166,7 @@ GResponse&amp; GResponse::operator=(const GResponse&amp; rsp)
  * @param[in] model Sky model.
  * @param[in] event Event.
  * @param[in] obs Observation.
- * @param[in] grad Should model gradients be computed? (default: true)
+ * @param[in] grad Should model gradients be computed?
  * @return Event probability.
  *
  * Computes the event probability
@@ -202,6 +203,9 @@ double GResponse::convolve(const GModelSky&amp;    model,
             // Retrieve true energy boundaries
             GEbounds ebounds = this->ebounds(event.energy());
  
+            // Initialise Ndarray array
+            GNdarray array;
+
             // Loop over all boundaries
             for (int i = 0; i < ebounds.size(); ++i) {
  
@@ -214,19 +218,67 @@ double GResponse::convolve(const GModelSky&amp;    model,
  
                     // Setup integration function
                     edisp_kern integrand(this, &obs, &model, &event, srcTime, grad);
-                    GIntegral  integral(&integrand);
+                    GIntegrals integral(&integrand);
  
                     // Set number of iterations
                     integral.fixed_iter(iter);
  
                     // Do Romberg integration
-                    prob += integral.romberg(etrue_min, etrue_max);
+                    if (array.size() == 0) {
+                        array = integral.romberg(etrue_min, etrue_max);
+                    }
+                    else {
+                        array += integral.romberg(etrue_min, etrue_max);
+                    }
  
                 } // endif: interval was valid
  
             } // endfor: looped over intervals
  
-        }
+            // Initialise array index
+            int index = 0;
+
+            // Get probability
+            prob = array(index++);
+
+            // Set gradients
+            if (grad) {
+
+                // Set spectral gradients
+                if (model.spectral() != NULL) {
+                    for (int i = 0; i < model.spectral()->size(); ++i) {
+                        GModelPar& par = (*(model.spectral()))[i];
+                        if (par.is_free() && par.has_grad()) {
+                            par.factor_gradient(array(index++));
+                        }
+                    }
+                }
+
+                // Set temporal gradients
+                if (model.temporal() != NULL) {
+                    for (int i = 0; i < model.temporal()->size(); ++i) {
+                        GModelPar& par = (*(model.temporal()))[i];
+                        if (par.is_free() && par.has_grad()) {
+                            par.factor_gradient(array(index++));
+                        }
+                    }
+                }
+
+                // Optionally set scale gradient for instrument
+                if (model.has_scales()) {
+                    for (int i = 0; i < model.scales(); ++i) {
+                        const GModelPar& par = model.scale(i);
+                        if (par.name() == obs.instrument()) {
+                            if (par.is_free() && par.has_grad()) {
+                                par.factor_gradient(array(index++));
+                            }
+                        }
+                    }
+                }
+
+            } // endif: set gradients
+
+        } // endif: Case A
  
         // Case B: No integration (assume no energy dispersion)
         else {
@@ -237,7 +289,7 @@ double GResponse::convolve(const GModelSky&amp;    model,
             // Evaluate probability
             prob = eval_prob(model, event, srcEng, srcTime, obs, grad);
  
-        }
+        } // endelse: Case B
  
         // Compile option: Check for NaN/Inf
         #if defined(G_NAN_CHECK)
@@ -1032,36 +1084,133 @@ double GResponse::eval_prob(const GModelSky&amp;    model,
  
  
 /***********************************************************************//**
- * @brief Integration kernel for GResponse::edisp_kern() class
+ * @brief Constructor for energy dispersion integration kernel class
+ *
+ * @param[in] parent Pointer to response.
+ * @param[in] obs Pointer to observation.
+ * @param[in] model Pointer to sky model.
+ * @param[in] event Pointer to event.
+ * @param[in] srcTime True time.
+ * @param[in] grad Evaluate gradients?
+ *
+ * This method constructs the integration kernel needed for the energy
+ * dispersion computation.
+ ***************************************************************************/
+GResponse::edisp_kern::edisp_kern(const GResponse*    parent,
+                                  const GObservation* obs,
+                                  const GModelSky*    model,
+                                  const GEvent*       event,
+                                  const GTime&        srcTime,
+                                  const bool&         grad)
+{
+    // Set members
+    m_parent  = parent;
+    m_obs     = obs;
+    m_model   = model;
+    m_event   = event;
+    m_srcTime = srcTime;
+    m_grad    = grad;
+
+    // If gradients are requested then put pointers to all relevant parameter
+    // in the parameter pointer vector
+    m_pars.clear();
+    if (m_grad) {
+        if (m_model->spectral() != NULL) {
+            for (int i = 0; i < m_model->spectral()->size(); ++i) {
+                GModelPar& par = (*(m_model->spectral()))[i];
+                if (par.is_free() && par.has_grad()) {
+                    m_pars.push_back(&par);
+                }
+            }
+        }
+        if (m_model->temporal() != NULL) {
+            for (int i = 0; i < m_model->temporal()->size(); ++i) {
+                GModelPar& par = (*(m_model->temporal()))[i];
+                if (par.is_free() && par.has_grad()) {
+                    m_pars.push_back(&par);
+                }
+            }
+        }
+        if (m_model->has_scales()) {
+            for (int i = 0; i < m_model->scales(); ++i) {
+                GModelPar& par = const_cast<GModelPar&>(m_model->scale(i));
+                if (par.name() == m_obs->instrument()) {
+                    if (par.is_free() && par.has_grad()) {
+                        m_pars.push_back(&par);
+                    }
+                }
+            }
+        }
+    }
+
+    // Allocate Ndarray to be returned as return type
+    m_array = GNdarray(m_pars.size()+1);
+
+    // Return
+    return;
+}
+
+
+/***********************************************************************//**
+ * @brief Return energy dispersion integration kernel
+ *
+ * @return Ndarray.
+ *
+ * This method returns an Ndarray with the last values that were computed
+ * by the eval() method. In case that eval() was not called before, the
+ * method returns an Ndarray with the right dimension and with all values
+ * set to zero.
+ ***************************************************************************/
+const GNdarray& GResponse::edisp_kern::array(void) const
+{
+    // Return Ndarray
+    return m_array;
+}
+
+
+/***********************************************************************//**
+ * @brief Evaluate energy dispersion integration kernel
  *
  * @param[in] etrue True photon energy in MeV.
  *
  * This method implements the integration kernel needed for the
  * GResponse::edisp_kern() class.
  ***************************************************************************/
-double GResponse::edisp_kern::eval(const double& etrue)
+GNdarray GResponse::edisp_kern::eval(const double& etrue)
 {
+    // Initialise Ndarray array index
+    int index = 0;
+
     // Set true energy
     GEnergy srcEng;
     srcEng.MeV(etrue);
  
     // Get function value
-    double value = m_parent->eval_prob(*m_model, *m_event,
-                                       srcEng, m_srcTime, *m_obs, m_grad);
+    m_array(index++) = m_parent->eval_prob(*m_model, *m_event,
+                                           srcEng, m_srcTime, *m_obs,
+                                           m_grad);
+
+    // If gradients are requested then extract them and put them into the
+    // array
+    if (m_grad) {
+        for (int i = 0; i < m_pars.size(); ++i, ++index) {
+            m_array(index) = m_pars[i]->factor_gradient();
+        }
+    }
  
     // Compile option: Check for NaN
     #if defined(G_NAN_CHECK)
-    if (gammalib::is_notanumber(value) || gammalib::is_infinite(value)) {
+    if (gammalib::is_notanumber(m_array(0)) || gammalib::is_infinite(m_array(0))) {
         std::cout << "*** ERROR: GResponse::edisp_kern::eval";
         std::cout << "(etrue=" << etrue << "): ";
         std::cout << " NaN/Inf encountered";
-        std::cout << " (value=" << value;
+        std::cout << " (value=" << m_array(0);
         std::cout << ")" << std::endl;
     }
     #endif
  
-    // Return value
-    return value;
+    // Return array
+    return m_array;
 }