From f54762c925814a56be31968153ee9ffa16d9f68b Mon Sep 17 00:00:00 2001 From: Jürgen Knödlseder Date: Sat, 8 Aug 2020 13:53:57 +0200 Subject: [PATCH] Implement specific cache for binned data at the level of GCTAResponseIrf (#3316) --- ChangeLog | 3 ++- NEWS | 8 ++++---- doc/source/admin/release_history/1.7.1.rst | 2 +- inst/cta/include/GCTAResponse.hpp | 10 ++++++++++ inst/cta/include/GCTAResponseCube.hpp | 10 ---------- inst/cta/pyext/GCTAResponse.i | 7 +++++-- inst/cta/pyext/GCTAResponseCube.i | 5 ----- inst/cta/pyext/GCTAResponseIrf.i | 4 +--- inst/cta/src/GCTAResponse.cpp | 203 +++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++ inst/cta/src/GCTAResponseCube.cpp | 192 ------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------ inst/cta/src/GCTAResponseIrf.cpp | 2 -- 11 files changed, 226 insertions(+), 220 deletions(-) diff --git a/ChangeLog b/ChangeLog index 1f6eb77..c34a6ed 100644 --- a/ChangeLog +++ b/ChangeLog @@ -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) diff --git a/NEWS b/NEWS index 7696c4a..d8907c0 100644 --- a/NEWS +++ b/NEWS @@ -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 diff --git a/doc/source/admin/release_history/1.7.1.rst b/doc/source/admin/release_history/1.7.1.rst index e9b159a..4e5f306 100644 --- a/doc/source/admin/release_history/1.7.1.rst +++ b/doc/source/admin/release_history/1.7.1.rst @@ -11,7 +11,7 @@ GammaLib 1.7.1 is a bug fix release for GammaLib release 1.7. The following bugs were fixed: * [`3316 `_] - - Implement specific response cache for CTA stacked analysis + Implement specific response cache for CTA binned and stacked analysis * [`3308 `_] - Speed up model definition XML file writing * [`3295 `_] - diff --git a/inst/cta/include/GCTAResponse.hpp b/inst/cta/include/GCTAResponse.hpp index ba27243..33055ba 100644 --- a/inst/cta/include/GCTAResponse.hpp +++ b/inst/cta/include/GCTAResponse.hpp @@ -29,6 +29,7 @@ /* __ Includes ___________________________________________________________ */ #include +#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 m_cache_names; //!< Model names + mutable std::vector m_cache_values; //!< Cached values }; #endif /* GCTARESPONSE_HPP */ diff --git a/inst/cta/include/GCTAResponseCube.hpp b/inst/cta/include/GCTAResponseCube.hpp index 8b9bb6b..21a0fe7 100644 --- a/inst/cta/include/GCTAResponseCube.hpp +++ b/inst/cta/include/GCTAResponseCube.hpp @@ -30,7 +30,6 @@ /* __ Includes ___________________________________________________________ */ #include #include -#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 m_cache_names; //!< Model names - mutable std::vector m_cache_values; //!< Cached values }; diff --git a/inst/cta/pyext/GCTAResponse.i b/inst/cta/pyext/GCTAResponse.i index ac0297d..b6d66ee 100644 --- a/inst/cta/pyext/GCTAResponse.i +++ b/inst/cta/pyext/GCTAResponse.i @@ -1,7 +1,7 @@ /*************************************************************************** * 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; }; diff --git a/inst/cta/pyext/GCTAResponseCube.i b/inst/cta/pyext/GCTAResponseCube.i index 4dd997d..c19015a 100644 --- a/inst/cta/pyext/GCTAResponseCube.i +++ b/inst/cta/pyext/GCTAResponseCube.i @@ -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); diff --git a/inst/cta/pyext/GCTAResponseIrf.i b/inst/cta/pyext/GCTAResponseIrf.i index b46f687..41bece6 100644 --- a/inst/cta/pyext/GCTAResponseIrf.i +++ b/inst/cta/pyext/GCTAResponseIrf.i @@ -1,7 +1,7 @@ /*************************************************************************** * 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; diff --git a/inst/cta/src/GCTAResponse.cpp b/inst/cta/src/GCTAResponse.cpp index b249780..69c40cb 100644 --- a/inst/cta/src/GCTAResponse.cpp +++ b/inst/cta/src/GCTAResponse.cpp @@ -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& GCTAResponse::operator=(const GCTAResponse& 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(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(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& GCTAResponse::operator=(const GCTAResponse& 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; } diff --git a/inst/cta/src/GCTAResponseCube.cpp b/inst/cta/src/GCTAResponseCube.cpp index d3d7cea..8184649 100644 --- a/inst/cta/src/GCTAResponseCube.cpp +++ b/inst/cta/src/GCTAResponseCube.cpp @@ -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& 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(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(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& 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; } diff --git a/inst/cta/src/GCTAResponseIrf.cpp b/inst/cta/src/GCTAResponseIrf.cpp index af61e4e..f0b0e47 100644 --- a/inst/cta/src/GCTAResponseIrf.cpp +++ b/inst/cta/src/GCTAResponseIrf.cpp @@ -2131,8 +2131,6 @@ void GCTAResponseIrf::copy_members(const GCTAResponseIrf& 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; -- libgit2 0.21.4