Implement specific response cache for CTA stacked analysis (#3316)

Jürgen Knödlseder
1 parent 30832424
Showing 7 changed files with 230 additions and 23 deletions
ChangeLog
NEWS
README.md
doc/source/admin/release_history/1.7.1.rst
inst/cta/include/GCTAResponseCube.hpp
inst/cta/pyext/GCTAResponseCube.i
inst/cta/src/GCTAResponseCube.cpp
-2020-08-07
+2020-08-08
  
         * Bug fix version 1.7.1 released
         ================================
  
+        Implement specific response cache for CTA stacked analysis (#3316)
         Speed up model definition XML file writing (#3308)
-        Switch off response cache for stacked CTA analysis (#3305)
         Fix automatic type casting in Python interface (#3295)
         Make parameter file reading and writing thread safe (#3287)
         Add quotes to sky region names (#3280)
 New Features and Important Changes in GammaLib 1.7.1
  
-7 August 2020
+8 August 2020
  
  
 1. Introduction
@@ -144,8 +144,10 @@ None
  
 20. CTA interface
 -----------------
-The response cache was switched off for CTA stacked analysis since it used up
-too much memory in case that many sources needed to be cached (#3305).
+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).
  
  
 21. COMPTEL interface
 GammaLib information
 ====================
-* Version: 1.7.1 (7 August 2020)
+* Version: 1.7.1 (8 August 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.7.1:
  
-GammaLib 1.7.1 release (7 August 2020)
+GammaLib 1.7.1 release (8 August 2020)
 ======================================
  
 Overview
@@ -10,10 +10,10 @@ 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
 * [`3308 <https://cta-redmine.irap.omp.eu/issues/3308>`_] -
   Speed up model definition XML file writing
-* [`3305 <https://cta-redmine.irap.omp.eu/issues/3305>`_] -
-  Switch off response cache for stacked CTA analysis
 * [`3295 <https://cta-redmine.irap.omp.eu/issues/3295>`_] -
   Fix automatic type casting in Python interface
 * [`3287 <https://cta-redmine.irap.omp.eu/issues/3287>`_] -
@@ -30,6 +30,7 @@
 /* __ Includes ___________________________________________________________ */
 #include <string>
 #include <vector>
+#include "GNdarray.hpp"
 #include "GCTAResponse.hpp"
 #include "GCTACubeExposure.hpp"
 #include "GCTACubePsf.hpp"
@@ -97,6 +98,11 @@ 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);
@@ -147,6 +153,10 @@ 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
 };
  
  
 /***************************************************************************
  *      GCTAResponseCube.i - CTA cube analysis response function class     *
  * ----------------------------------------------------------------------- *
- *  copyright (C) 2014-2018 by Juergen Knoedlseder                         *
+ *  copyright (C) 2014-2020 by Juergen Knoedlseder                         *
  * ----------------------------------------------------------------------- *
  *                                                                         *
  *  This program is free software: you can redistribute it and/or modify   *
@@ -68,6 +68,11 @@ 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);
@@ -31,17 +31,6 @@
 #include <cmath>
 #include <string>
 #include "GTools.hpp"
-#include "GCTAResponseCube.hpp"
-#include "GCTAResponse_helpers.hpp"
-#include "GCTAInstDir.hpp"
-#include "GCTAEventBin.hpp"
-#include "GCTASupport.hpp"
-#include "GModelSpatialPointSource.hpp"
-#include "GModelSpatialRadial.hpp"
-#include "GModelSpatialRadialShell.hpp"
-#include "GModelSpatialElliptical.hpp"
-#include "GModelSpatialComposite.hpp"
-#include "GModelSpatialDiffuse.hpp"
 #include "GPhoton.hpp"
 #include "GSource.hpp"
 #include "GEvent.hpp"
@@ -50,6 +39,19 @@
 #include "GTime.hpp"
 #include "GIntegral.hpp"
 #include "GObservation.hpp"
+#include "GModelSpatialPointSource.hpp"
+#include "GModelSpatialRadial.hpp"
+#include "GModelSpatialRadialShell.hpp"
+#include "GModelSpatialElliptical.hpp"
+#include "GModelSpatialComposite.hpp"
+#include "GModelSpatialDiffuse.hpp"
+#include "GCTAResponseCube.hpp"
+#include "GCTAResponse_helpers.hpp"
+#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&)"
@@ -443,6 +445,188 @@ 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.
@@ -635,8 +819,10 @@ void GCTAResponseCube::init_members(void)
     m_background.clear();
     m_apply_edisp    = false;
     m_has_edisp      = false;
-    m_use_irf_cache  = false;  //!< Current response cache uses too
-    m_use_nroi_cache = false;  //!< much memory (#3305)
+
+    // Kludge: Initialise cube response cache
+    m_cache_names.clear();
+    m_cache_values.clear();
  
     // Return
     return;
@@ -658,6 +844,10 @@ 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;
 }