Implement specific cache for binned data at the level of GCTAResponseIrf (#3316)

By moving the specific cache from GCTAReponseCube to GResponseIrf, the specific cache is also available for binned analysis.

Implement specific cache for binned data at the level of GCTAResponseIrf (#3316)
By moving the specific cache from GCTAReponseCube to GResponseIrf, the specific cache is also available for binned analysis.
Jürgen Knödlseder
1 parent bc5cda86
Showing 11 changed files with 226 additions and 220 deletions
ChangeLog
NEWS
doc/source/admin/release_history/1.7.1.rst
inst/cta/include/GCTAResponse.hpp
inst/cta/include/GCTAResponseCube.hpp
inst/cta/pyext/GCTAResponse.i
inst/cta/pyext/GCTAResponseCube.i
inst/cta/pyext/GCTAResponseIrf.i
inst/cta/src/GCTAResponse.cpp
inst/cta/src/GCTAResponseCube.cpp
inst/cta/src/GCTAResponseIrf.cpp
@@ -3,7 +3,8 @@
         * Bug fix version 1.7.1 released
         ================================
-        Implement specific response cache for CTA stacked analysis (#3316)
+        Implement specific response cache for CTA binned and stacked
+          analysis (#3316)
         Speed up model definition XML file writing (#3308)
         Fix automatic type casting in Python interface (#3295)
         Make parameter file reading and writing thread safe (#3287)
@@ -144,10 +144,10 @@ None
 20. CTA interface
 -----------------
-A specific response cache was implemented for the stacked analysis that reduces
-the memory requirements by at least a factor of 6. Instead of using the std::map
-cache, the full data-space is cached, avoiding the need for any bookkeeping
-(#3316).
+A specific response cache was implemented for the binned and stacked analysis
+that reduces the memory requirements by at least a factor of 6. Instead of
+using the std::map cache, the full data cube is cached, avoiding the need for
+any bookkeeping that takes a lot of memory in the traditional cache (#3316).
 21. COMPTEL interface
@@ -11,7 +11,7 @@ GammaLib 1.7.1 is a bug fix release for GammaLib release 1.7.
 The following bugs were fixed:
 * [`3316 <https://cta-redmine.irap.omp.eu/issues/3316>`_] -
-  Implement specific response cache for CTA stacked analysis
+  Implement specific response cache for CTA binned and stacked analysis
 * [`3308 <https://cta-redmine.irap.omp.eu/issues/3308>`_] -
   Speed up model definition XML file writing
 * [`3295 <https://cta-redmine.irap.omp.eu/issues/3295>`_] -
@@ -29,6 +29,7 @@
 /* __ Includes ___________________________________________________________ */
 #include <string>
+#include "GNdarray.hpp"
 #include "GResponse.hpp"
 /* __ Type definitions ___________________________________________________ */
@@ -84,11 +85,20 @@ public:
     virtual void          write(GXmlElement& xml) const = 0;
     virtual std::string   print(const GChatter& chatter = NORMAL) const = 0;
+    // Overloaded virtual base class methods
+    virtual double irf_spatial(const GEvent&       event,
+                               const GSource&      source,
+                               const GObservation& obs) const;
+
 protected:
     // Protected methods
     void init_members(void);
     void copy_members(const GCTAResponse& rsp);
     void free_members(void);
+
+    // Kludge: cube response cache
+    mutable std::vector<std::string> m_cache_names;  //!< Model names
+    mutable std::vector<GNdarray>    m_cache_values; //!< Cached values
 };
 #endif /* GCTARESPONSE_HPP */
@@ -30,7 +30,6 @@
 /* __ Includes ___________________________________________________________ */
 #include <string>
 #include <vector>
-#include "GNdarray.hpp"
 #include "GCTAResponse.hpp"
 #include "GCTACubeExposure.hpp"
 #include "GCTACubePsf.hpp"
@@ -98,11 +97,6 @@ public:
     virtual void              write(GXmlElement& xml) const;
     virtual std::string       print(const GChatter& chatter = NORMAL) const;
-    // Overloaded virtual base class methods
-    virtual double      irf_spatial(const GEvent&       event,
-                                    const GSource&      source,
-                                    const GObservation& obs) const;
-
     // Other Methods
     const GCTACubeExposure&   exposure(void) const;
     void                      exposure(const GCTACubeExposure& exposure);
@@ -153,10 +147,6 @@ private:
     GCTACubeEdisp      m_edisp;       //!< Energy dispersion cube
     mutable bool       m_apply_edisp; //!< Apply energy dispersion
     mutable bool       m_has_edisp;   //!< Flag to indicate if energy
-
-    // Kludge: cube response cache
-    mutable std::vector<std::string> m_cache_names;  //!< Model names
-    mutable std::vector<GNdarray>    m_cache_values; //!< Cached values
 };
 /***************************************************************************
  *             GCTAResponse.i - CTA response abstract base class           *
  * ----------------------------------------------------------------------- *
- *  copyright (C) 2010-2019 by Juergen Knoedlseder                         *
+ *  copyright (C) 2010-2020 by Juergen Knoedlseder                         *
  * ----------------------------------------------------------------------- *
  *                                                                         *
  *  This program is free software: you can redistribute it and/or modify   *
@@ -58,7 +58,10 @@ public:
     virtual void          read(const GXmlElement& xml) = 0;
     virtual void          write(GXmlElement& xml) const = 0;
-    // Implemeted methods
+    // Overloaded virtual base class methods
+    virtual double irf_spatial(const GEvent&       event,
+                               const GSource&      source,
+                               const GObservation& obs) const;
 };
@@ -68,11 +68,6 @@ public:
     virtual void              read(const GXmlElement& xml);
     virtual void              write(GXmlElement& xml) const;
-    // Overloaded virtual base class methods
-    virtual double      irf_spatial(const GEvent&       event,
-                                    const GSource&      source,
-                                    const GObservation& obs) const;
-
     // Other Methods
     const GCTACubeExposure&   exposure(void) const;
     void                      exposure(const GCTACubeExposure& exposure);
 /***************************************************************************
  *        GCTAResponseIrf.i - CTA instrument response function class       *
  * ----------------------------------------------------------------------- *
- *  copyright (C) 2010-2018 by Juergen Knoedlseder                         *
+ *  copyright (C) 2010-2020 by Juergen Knoedlseder                         *
  * ----------------------------------------------------------------------- *
  *                                                                         *
  *  This program is free software: you can redistribute it and/or modify   *
@@ -60,8 +60,6 @@ public:
     virtual void             read(const GXmlElement& xml);
     virtual void             write(GXmlElement& xml) const;
-    // Overload virtual base class methods
-
     // Other Methods
     bool                  apply_edisp(void) const;
     void                  apply_edisp(const bool& apply_edisp) const;
@@ -33,6 +33,8 @@
 #include "GObservations.hpp"
 #include "GModelSpatial.hpp"
 #include "GModelSpatialComposite.hpp"
+#include "GCTAEventCube.hpp"               // Kludge
+#include "GCTAEventBin.hpp"                // Kludge
 #include "GCTAResponse.hpp"
 #include "GCTASupport.hpp"
@@ -151,6 +153,199 @@ GCTAResponse&amp; GCTAResponse::operator=(const GCTAResponse&amp; rsp)
  =                                                                         =
  ==========================================================================*/
+/***********************************************************************//**
+ * @brief Return instrument response integrated over the spatial model
+ *
+ * @param[in] event Event.
+ * @param[in] source Source.
+ * @param[in] obs Observation.
+ * @return Instrument response to a spatial model.
+ *
+ * Returns the instrument response for a given event, source and observation
+ * integrated over the spatial model component. The method computes
+ *
+ * \f[
+ *    {\tt irf}(p', E', t') = \int_p M_{\rm S}(p | E, t) \,
+ *                                   R(p', E', t' | p, E, t) \, d\,p
+ * \f]
+ *
+ * where
+ * * \f$M_{\rm S}(p | E, t)\f$ is the spatial model component,
+ * * \f$R(p', E', t' | p, E, t)\f$ is the Instrument Response Function (IRF),
+ * * \f$p'\f$ is the measured instrument direction,
+ * * \f$E'\f$ is the measured or reconstructed energy,
+ * * \f$t'\f$ is the measured arrival time,
+ * * \f$p\f$ is the true photon arrival direction,
+ * * \f$E\f$ is the true photon energy, and
+ * * \f$t\f$ is the true trigger time.
+ *
+ * The integration is done over all relevant true sky directions \f$p\f$.
+ *
+ * Depending on the type of the source model the method branches to the
+ * following methods to perform the actual computations
+ *
+ *      irf_ptsrc() - for the handling of a point source
+ *      irf_radial() - for radial models
+ *      irf_elliptical() - for elliptical models
+ *      irf_diffuse() - for diffuse models
+ *      irf_composite() - for composite models
+ *
+ * The method implements a caching mechanism for spatial models that have all
+ * parameters fixed. For those models the instrument response for a given
+ * event and observation is only computed once and then stored in an internal
+ * cache from which it is fetched back in case that the method is called
+ * again for the same event and observation.
+ ***************************************************************************/
+double GCTAResponse::irf_spatial(const GEvent&       event,
+                                 const GSource&      source,
+                                 const GObservation& obs) const
+{
+    // Initialise IRF value
+    double irf = 0.0;
+
+    // Set IRF value attributes
+    std::string     name  = obs.id() + "::" + source.name();
+    const GInstDir& dir   = event.dir();
+    const GEnergy&  ereco = event.energy();
+    const GEnergy&  etrue = source.energy();
+
+    // Signal if spatial model has free parameters
+    bool has_free_pars = source.model()->has_free_pars();
+
+    // Kludge: if the response cache should be used, the event is a bin, the
+    // model has no free parameters, and there is no energy dispersion then
+    // use the special cache for binned or stacked analysis
+    bool use_kludge = false;
+    if (m_use_irf_cache && event.is_bin() && !has_free_pars && !use_edisp()) {
+
+        // Get reference to CTA event bin
+        const GCTAEventBin& bin = static_cast<const GCTAEventBin&>(event);
+
+        // If pixel and energy indices are valid then proceed with the kludge
+        if (bin.ipix() >= 0 && bin.ieng() >=0) {
+
+            // Signal that kludge will be used
+            use_kludge = true;
+
+            // Build unique cache name
+            std::string name  = obs.id() + "::" + source.name();
+
+            // Get index in cache, returns -1 if name is not found in cache
+            int index = -1;
+            for (int i = 0; i < m_cache_names.size(); ++i) {
+                if (m_cache_names[i] == name) {
+                    index = i;
+                    break;
+                }
+            }
+
+            // If index was not found then allocate a new cache map
+            if (index == -1) {
+
+                // Get pointer to event cube
+                const GCTAEventCube* cube =
+                      static_cast<const GCTAEventCube*>(obs.events());
+
+                // Allocate cache
+                GNdarray cache(cube->nx()*cube->ny(), cube->ebins());
+
+                // Initialise all cache values with -1 (not set)
+                for (int i = 0; i < cache.size(); ++i) {
+                    cache(i) = -1.0;
+                }
+
+                // Insert cache
+                m_cache_names.push_back(name);
+                m_cache_values.push_back(cache);
+
+                // Set index
+                index = m_cache_names.size()-1;
+
+            } // endif: allocated new cache
+
+            // Get cache value
+            irf = m_cache_values[index](bin.ipix(), bin.ieng());
+
+            // If cache value is not valid then compute IRF
+            if (irf < 0.0) {
+
+                // Compute IRF for spatial model
+                switch (source.model()->code()) {
+                    case GMODEL_SPATIAL_POINT_SOURCE:
+                        irf = irf_ptsrc(event, source, obs);
+                        break;
+                    case GMODEL_SPATIAL_RADIAL:
+                        irf = irf_radial(event, source, obs);
+                        break;
+                    case GMODEL_SPATIAL_ELLIPTICAL:
+                        irf = irf_elliptical(event, source, obs);
+                        break;
+                    case GMODEL_SPATIAL_DIFFUSE:
+                        irf = irf_diffuse(event, source, obs);
+                        break;
+                    case GMODEL_SPATIAL_COMPOSITE:
+                        irf = irf_composite(event, source, obs);
+                        break;
+                    default:
+                        break;
+                }
+
+                // Set cache value
+                m_cache_values[index](bin.ipix(), bin.ieng()) = irf;
+
+            } // endif: computed IRF
+
+        } // endif: event bin was part of cube
+
+    } // endif: kludge
+
+    // If kludge was not used then use the release 1.7 response cache
+    if (!use_kludge) {
+
+        // If the spatial model component has free parameters, or the response
+        // cache should not be used, or the cache does not contain the requested
+        // IRF value then compute the IRF value for the spatial model.
+        if (has_free_pars    ||
+            !m_use_irf_cache ||
+            !m_irf_cache.contains(name, dir, ereco, etrue, &irf)) {
+
+            // Compute IRF for spatial model
+            switch (source.model()->code()) {
+                case GMODEL_SPATIAL_POINT_SOURCE:
+                    irf = irf_ptsrc(event, source, obs);
+                    break;
+                case GMODEL_SPATIAL_RADIAL:
+                    irf = irf_radial(event, source, obs);
+                    break;
+                case GMODEL_SPATIAL_ELLIPTICAL:
+                    irf = irf_elliptical(event, source, obs);
+                    break;
+                case GMODEL_SPATIAL_DIFFUSE:
+                    irf = irf_diffuse(event, source, obs);
+                    break;
+                case GMODEL_SPATIAL_COMPOSITE:
+                    irf = irf_composite(event, source, obs);
+                    break;
+                default:
+                    break;
+            }
+
+        } // endif: computed spatial model
+
+        // If the spatial model has no free parameters and the response cache
+        // should be used then put the IRF value in the response cache.
+        if (!has_free_pars && m_use_irf_cache) {
+            m_irf_cache.set(name, dir, ereco, etrue, irf);
+        }
+
+    } // endelse: used release 1.7 response cache
+
+    // Return IRF value
+    return irf;
+}
+
+
+
 /*==========================================================================
  =                                                                         =
  =                             Private methods                             =
@@ -162,6 +357,10 @@ GCTAResponse&amp; GCTAResponse::operator=(const GCTAResponse&amp; rsp)
  ***************************************************************************/
 void GCTAResponse::init_members(void)
 {
+    // Kludge: Initialise cube response cache
+    m_cache_names.clear();
+    m_cache_values.clear();
+
     // Return
     return;
 }
@@ -174,6 +373,10 @@ void GCTAResponse::init_members(void)
  ***************************************************************************/
 void GCTAResponse::copy_members(const GCTAResponse& rsp)
 {
+    // Kludge: Copy cube response cache
+    m_cache_names  = rsp.m_cache_names;
+    m_cache_values = rsp.m_cache_values;
+
     // Return
     return;
 }
@@ -50,8 +50,6 @@
 #include "GCTAInstDir.hpp"
 #include "GCTAEventBin.hpp"
 #include "GCTASupport.hpp"
-#include "GCTAEventCube.hpp"               // Kludge
-#include "GCTAEventBin.hpp"                // Kludge
 /* __ Method name definitions ____________________________________________ */
 #define G_IRF        "GCTAResponseCube::irf(GEvent&, GPhoton& GObservation&)"
@@ -445,188 +443,6 @@ GEbounds GCTAResponseCube::ebounds(const GEnergy&amp; obsEng) const
 /***********************************************************************//**
- * @brief Return instrument response integrated over the spatial model
- *
- * @param[in] event Event.
- * @param[in] source Source.
- * @param[in] obs Observation.
- * @return Instrument response to a spatial model.
- *
- * Returns the instrument response for a given event, source and observation
- * integrated over the spatial model component. The method computes
- *
- * \f[
- *    {\tt irf}(p', E', t') = \int_p M_{\rm S}(p | E, t) \,
- *                                   R(p', E', t' | p, E, t) \, d\,p
- * \f]
- *
- * where
- * * \f$M_{\rm S}(p | E, t)\f$ is the spatial model component,
- * * \f$R(p', E', t' | p, E, t)\f$ is the Instrument Response Function (IRF),
- * * \f$p'\f$ is the measured instrument direction,
- * * \f$E'\f$ is the measured or reconstructed energy,
- * * \f$t'\f$ is the measured arrival time,
- * * \f$p\f$ is the true photon arrival direction,
- * * \f$E\f$ is the true photon energy, and
- * * \f$t\f$ is the true trigger time.
- *
- * The integration is done over all relevant true sky directions \f$p\f$.
- *
- * Depending on the type of the source model the method branches to the
- * following methods to perform the actual computations
- *
- *      irf_ptsrc() - for the handling of a point source
- *      irf_radial() - for radial models
- *      irf_elliptical() - for elliptical models
- *      irf_diffuse() - for diffuse models
- *      irf_composite() - for composite models
- *
- * The method implements a caching mechanism for spatial models that have all
- * parameters fixed. For those models the instrument response for a given
- * event and observation is only computed once and then stored in an internal
- * cache from which it is fetched back in case that the method is called
- * again for the same event and observation.
- ***************************************************************************/
-double GCTAResponseCube::irf_spatial(const GEvent&       event,
-                                     const GSource&      source,
-                                     const GObservation& obs) const
-{
-    // Initialise IRF value
-    double irf = 0.0;
-
-    // Set IRF value attributes
-    std::string     name  = obs.id() + "::" + source.name();
-    const GInstDir& dir   = event.dir();
-    const GEnergy&  ereco = event.energy();
-    const GEnergy&  etrue = source.energy();
-
-    // Signal if spatial model has free parameters
-    bool has_free_pars = source.model()->has_free_pars();
-
-    // Kludge: the IRF response cache should be used, the model has no
-    // free parameters and there is no energy dispersion
-    if (m_use_irf_cache && !has_free_pars && !use_edisp()) {
-
-        // Build unique cache name
-        std::string name  = obs.id() + "::" + source.name();
-
-        // Get index in cache, returns -1 if name is not found in cache
-        int index = -1;
-        for (int i = 0; i < m_cache_names.size(); ++i) {
-            if (m_cache_names[i] == name) {
-                index = i;
-                break;
-            }
-        }
-
-        // If index was not found then allocate a new cache map
-        if (index == -1) {
-
-            // Get pointer to event cube
-            const GCTAEventCube* cube =
-                  static_cast<const GCTAEventCube*>(obs.events());
-
-            // Allocate cache
-            GNdarray cache(cube->nx()*cube->ny(), cube->ebins());
-
-            // Initialise all cache values with -1 (not set)
-            for (int i = 0; i < cache.size(); ++i) {
-                cache(i) = -1.0;
-            }
-
-            // Insert cache
-            m_cache_names.push_back(name);
-            m_cache_values.push_back(cache);
-
-            // Set index
-            index = m_cache_names.size()-1;
-
-        } // endif: allocated new cache
-
-        // Get reference to CTA event bin
-        const GCTAEventBin& bin = static_cast<const GCTAEventBin&>(event);
-
-        // Get cache value
-        irf = m_cache_values[index](bin.ipix(), bin.ieng());
-
-        // If cache value is not valid then copute IRF
-        if (irf < 0.0) {
-
-            // Compute IRF for spatial model
-            switch (source.model()->code()) {
-                case GMODEL_SPATIAL_POINT_SOURCE:
-                    irf = irf_ptsrc(event, source, obs);
-                    break;
-                case GMODEL_SPATIAL_RADIAL:
-                    irf = irf_radial(event, source, obs);
-                    break;
-                case GMODEL_SPATIAL_ELLIPTICAL:
-                    irf = irf_elliptical(event, source, obs);
-                    break;
-                case GMODEL_SPATIAL_DIFFUSE:
-                    irf = irf_diffuse(event, source, obs);
-                    break;
-                case GMODEL_SPATIAL_COMPOSITE:
-                    irf = irf_composite(event, source, obs);
-                    break;
-                default:
-                    break;
-            }
-
-            // Set cache value
-            m_cache_values[index](bin.ipix(), bin.ieng()) = irf;
-
-        } // endif: computed IRF
-
-    } // endif: kludge
-
-    // ... otherwise use release 1.7 response cache
-    else {
-
-        // If the spatial model component has free parameters, or the response
-        // cache should not be used, or the cache does not contain the requested
-        // IRF value then compute the IRF value for the spatial model.
-        if (has_free_pars    ||
-            !m_use_irf_cache ||
-            !m_irf_cache.contains(name, dir, ereco, etrue, &irf)) {
-
-            // Compute IRF for spatial model
-            switch (source.model()->code()) {
-                case GMODEL_SPATIAL_POINT_SOURCE:
-                    irf = irf_ptsrc(event, source, obs);
-                    break;
-                case GMODEL_SPATIAL_RADIAL:
-                    irf = irf_radial(event, source, obs);
-                    break;
-                case GMODEL_SPATIAL_ELLIPTICAL:
-                    irf = irf_elliptical(event, source, obs);
-                    break;
-                case GMODEL_SPATIAL_DIFFUSE:
-                    irf = irf_diffuse(event, source, obs);
-                    break;
-                case GMODEL_SPATIAL_COMPOSITE:
-                    irf = irf_composite(event, source, obs);
-                    break;
-                default:
-                    break;
-            }
-
-        } // endif: computed spatial model
-
-        // If the spatial model has no free parameters and the response cache
-        // should be used then put the IRF value in the response cache.
-        if (!has_free_pars && m_use_irf_cache) {
-            m_irf_cache.set(name, dir, ereco, etrue, irf);
-        }
-
-    } // endelse: used release 1.7 response cache
-
-    // Return IRF value
-    return irf;
-}
-
-
-/***********************************************************************//**
  * @brief Read response information from XML element
  *
  * @param[in] xml XML element.
@@ -820,10 +636,6 @@ void GCTAResponseCube::init_members(void)
     m_apply_edisp    = false;
     m_has_edisp      = false;
-    // Kludge: Initialise cube response cache
-    m_cache_names.clear();
-    m_cache_values.clear();
-
     // Return
     return;
 }
@@ -844,10 +656,6 @@ void GCTAResponseCube::copy_members(const GCTAResponseCube&amp; rsp)
     m_apply_edisp = rsp.m_apply_edisp;
     m_has_edisp   = rsp.m_has_edisp;
-    // Kludge: Copy cube response cache
-    m_cache_names  = rsp.m_cache_names;
-    m_cache_values = rsp.m_cache_values;
-
     // Return
     return;
 }
@@ -2131,8 +2131,6 @@ void GCTAResponseIrf::copy_members(const GCTAResponseIrf&amp; rsp)
     m_xml_caldb      = rsp.m_xml_caldb;
     m_xml_rspname    = rsp.m_xml_rspname;
-    // Copy nroi cache
-
     // Clone members
     m_aeff       = (rsp.m_aeff       != NULL) ? rsp.m_aeff->clone()  : NULL;
     m_psf        = (rsp.m_psf        != NULL) ? rsp.m_psf->clone()   : NULL;