From d3768de7b6230a91048b26ce3bb52a6164d216e6 Mon Sep 17 00:00:00 2001 From: Jürgen Knödlseder Date: Tue, 20 Oct 2020 16:11:32 +0200 Subject: [PATCH] Implement sparse vector response cache (#3403) --- ChangeLog | 3 ++- NEWS | 7 ++++++- README.md | 2 +- doc/source/admin/release_history/2.0.rst | 4 +++- include/GResponse.hpp | 6 ++++-- include/GResponseVectorCache.hpp | 130 ++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++ include/GammaLib.hpp | 1 + include/Makefile.am | 1 + pyext/GResponseVectorCache.i | 78 ++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++ pyext/gammalib/obs.i | 1 + src/obs/GResponse.cpp | 66 ++++++++++++++++++++++++++++++++++++++++++++---------------------- src/obs/GResponseVectorCache.cpp | 480 ++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++ src/obs/Makefile.am | 1 + test/test_GObservation.cpp | 84 ++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++-- test/test_GObservation.hpp | 3 ++- 15 files changed, 836 insertions(+), 31 deletions(-) create mode 100644 include/GResponseVectorCache.hpp create mode 100644 pyext/GResponseVectorCache.i create mode 100644 src/obs/GResponseVectorCache.cpp diff --git a/ChangeLog b/ChangeLog index 9e5b44a..bd784c7 100644 --- a/ChangeLog +++ b/ChangeLog @@ -1,8 +1,9 @@ -2020-10-17 +2020-10-20 * Version 2.0.0 released ======================== + Implement sparse vector response cache (#3403) Implement GMatrixSparse::row() and GMatrixSparse::add_to_row() methods Integrate analytical gradients for energy dispersion computation (#3333) Add GFunctions abstract base class (#3328) diff --git a/NEWS b/NEWS index 107d1ee..e2b5a5f 100644 --- a/NEWS +++ b/NEWS @@ -1,6 +1,6 @@ New Features and Important Changes in GammaLib 2.0.0 -17 October 2020 +20 October 2020 1. Introduction @@ -14,6 +14,7 @@ made since the last release of GammaLib. The following classes have been added: - GFunctions - GIntegrals +- GResponseVectorCache The following classes have been removed: - none @@ -87,6 +88,10 @@ of functions (#3328). 10. Observation module ---------------------- +The GResponseVectorCache class was added to handle caching of the vector response. +Only non-zero IRF values are stored by the class. The class was integrated in the +GResponse class (#3403). + A vectorised response computation was implemented in GResponse that computes the IRF for all events in a given observations. This allows for a considerable reduction of the operations, in particular for binned or stacked analysis, where the code can diff --git a/README.md b/README.md index c2c1dd3..3e31858 100644 --- a/README.md +++ b/README.md @@ -1,6 +1,6 @@ GammaLib information ==================== -* Version: 2.0.0.dev (17 October 2020) +* Version: 2.0.0.dev (20 October 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/) diff --git a/doc/source/admin/release_history/2.0.rst b/doc/source/admin/release_history/2.0.rst index 13a1e1b..8a26fe7 100644 --- a/doc/source/admin/release_history/2.0.rst +++ b/doc/source/admin/release_history/2.0.rst @@ -1,6 +1,6 @@ .. _2.0: -GammaLib 2.0 release (17 October 2020) +GammaLib 2.0 release (20 October 2020) ====================================== Overview @@ -31,6 +31,8 @@ Improvements New features ------------ +* [`3403 `_] - + Implement sparse vector response cache * [`3328 `_] - Add ``GFunctions`` abstract base class and ``GIntegrals`` class diff --git a/include/GResponse.hpp b/include/GResponse.hpp index c1f28df..54fef63 100644 --- a/include/GResponse.hpp +++ b/include/GResponse.hpp @@ -33,6 +33,7 @@ #include "GFunction.hpp" #include "GTime.hpp" #include "GResponseCache.hpp" +#include "GResponseVectorCache.hpp" #include "GFunctions.hpp" /* __ Forward declarations _______________________________________________ */ @@ -319,8 +320,9 @@ protected: // Cache for irf(GEvent&, GSource&, GObservation&) and // nroi(GModelSky&, GEnergy&, GTime&, GEnergy&, GTime&, GObservation&) // computations - mutable GResponseCache m_irf_cache; - mutable GResponseCache m_nroi_cache; + mutable GResponseCache m_irf_cache; + mutable GResponseCache m_nroi_cache; + mutable GResponseVectorCache m_irf_vector_cache; }; #endif /* GRESPONSE_HPP */ diff --git a/include/GResponseVectorCache.hpp b/include/GResponseVectorCache.hpp new file mode 100644 index 0000000..87a25f6 --- /dev/null +++ b/include/GResponseVectorCache.hpp @@ -0,0 +1,130 @@ +/*************************************************************************** + * GResponseVectorCache.hpp - Response vector cache class * + * ----------------------------------------------------------------------- * + * copyright (C) 2020 by Juergen Knoedlseder * + * ----------------------------------------------------------------------- * + * * + * This program is free software: you can redistribute it and/or modify * + * it under the terms of the GNU General Public License as published by * + * the Free Software Foundation, either version 3 of the License, or * + * (at your option) any later version. * + * * + * This program is distributed in the hope that it will be useful, * + * but WITHOUT ANY WARRANTY; without even the implied warranty of * + * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the * + * GNU General Public License for more details. * + * * + * You should have received a copy of the GNU General Public License * + * along with this program. If not, see . * + * * + ***************************************************************************/ +/** + * @file GResponseVectorCache.hpp + * @brief Response vector cache class definition + * @author Juergen Knoedlseder + */ + +#ifndef GRESPONSEVECTORCACHE_HPP +#define GRESPONSEVECTORCACHE_HPP + +/* __ Includes ___________________________________________________________ */ +#include +#include +#include "GBase.hpp" + +/* __ Forward declarations _______________________________________________ */ +class GVector; + +/* __ Constants __________________________________________________________ */ + +/* __ Typedefs ___________________________________________________________ */ + + +/***********************************************************************//** + * @class GResponseVectorCache + * + * @brief Response vector cache class + * + * The class implements a vector cache for the Instrument Response Function + * values so that the values do not need to be recomputed each time but can + * be fetched from the cache. + ***************************************************************************/ +class GResponseVectorCache : public GBase { + +public: + // Constructors and destructors + GResponseVectorCache(void); + GResponseVectorCache(const GResponseVectorCache& cache); + virtual ~GResponseVectorCache(void); + + // Operators + GResponseVectorCache& operator=(const GResponseVectorCache& cache); + + // Methods + void clear(void); + GResponseVectorCache* clone(void) const; + std::string classname(void) const; + bool is_empty(void) const; + int size(void) const; + void set(const std::string& cache_id, + const GVector& vector); + void remove(const std::string& cache_id); + bool contains(const std::string& cache_id, + GVector* irfs = NULL) const; + std::string print(const GChatter& chatter = NORMAL) const; + +protected: + // Protected methods + void init_members(void); + void copy_members(const GResponseVectorCache& cache); + void free_members(void); + int find_cache(const std::string& cache_id) const; + + // Protected members + std::vector m_cache_ids; + std::vector m_cache_entries; + std::vector m_cache_values; + std::vector m_cache_indices; +}; + + +/***********************************************************************//** + * @brief Return class name + * + * @return String containing the class name ("GResponseVectorCache"). + ***************************************************************************/ +inline +std::string GResponseVectorCache::classname(void) const +{ + return ("GResponseVectorCache"); +} + + +/***********************************************************************//** + * @brief Checks whether the cache is empty + * + * @return True if cache is empty, false otherwise. + * + * Checks whether the response cache is empty. + ***************************************************************************/ +inline +bool GResponseVectorCache::is_empty(void) const +{ + return (m_cache_ids.empty()); +} + + +/***********************************************************************//** + * @brief Remove cache for source + * + * @param[in] name Source name. + * + * Remove cache for source with @p name. + ***************************************************************************/ +inline +int GResponseVectorCache::size(void) const +{ + return (m_cache_ids.size()); +} + +#endif /* GRESPONSEVECTORCACHE_HPP */ diff --git a/include/GammaLib.hpp b/include/GammaLib.hpp index f4f5895..e5f7a41 100644 --- a/include/GammaLib.hpp +++ b/include/GammaLib.hpp @@ -191,6 +191,7 @@ #include "GInstDir.hpp" #include "GResponse.hpp" #include "GResponseCache.hpp" +#include "GResponseVectorCache.hpp" #include "GPhotons.hpp" #include "GPhoton.hpp" #include "GSource.hpp" diff --git a/include/Makefile.am b/include/Makefile.am index 734835a..2a8ec44 100644 --- a/include/Makefile.am +++ b/include/Makefile.am @@ -151,6 +151,7 @@ pkginclude_HEADERS = GammaLib.hpp \ GInstDir.hpp \ GResponse.hpp \ GResponseCache.hpp \ + GResponseVectorCache.hpp \ GPhotons.hpp \ GPhoton.hpp \ GSource.hpp \ diff --git a/pyext/GResponseVectorCache.i b/pyext/GResponseVectorCache.i new file mode 100644 index 0000000..7aaf501 --- /dev/null +++ b/pyext/GResponseVectorCache.i @@ -0,0 +1,78 @@ +/*************************************************************************** + * GResponseVectorCache.i - Response vector cache class * + * ----------------------------------------------------------------------- * + * copyright (C) 2020 by Juergen Knoedlseder * + * ----------------------------------------------------------------------- * + * * + * This program is free software: you can redistribute it and/or modify * + * it under the terms of the GNU General Public License as published by * + * the Free Software Foundation, either version 3 of the License, or * + * (at your option) any later version. * + * * + * This program is distributed in the hope that it will be useful, * + * but WITHOUT ANY WARRANTY; without even the implied warranty of * + * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the * + * GNU General Public License for more details. * + * * + * You should have received a copy of the GNU General Public License * + * along with this program. If not, see . * + * * + ***************************************************************************/ +/** + * @file GResponseVectorCache.i + * @brief Response vector cache class definition + * @author Juergen Knoedlseder + */ +%{ +/* Put headers and other declarations here that are needed for compilation */ +#include "GResponseVectorCache.hpp" +%} + +/* __ Typemaps ___________________________________________________________ */ +/* + * Typemap so that the contains() method returns two values, first the + * flag indicating whether the value exists or not, and second the + * value, i.e. + * + * flag, value = cache.contains('Crab', GEnergy(1.0,'TeV'), GEnergy(1.0,'TeV')) + * + */ +%include "typemaps.i" +%apply GVector *OUTPUT { GVector* value }; + + +/***********************************************************************//** + * @class GResponseVectorCache + * + * @brief Response cache class + ***************************************************************************/ +class GResponseVectorCache : public GBase { + +public: + // Constructors and destructors + GResponseVectorCache(void); + GResponseVectorCache(const GResponseVectorCache& cache); + virtual ~GResponseVectorCache(void); + + // Methods + void clear(void); + GResponseVectorCache* clone(void) const; + std::string classname(void) const; + bool is_empty(void) const; + int size(void) const; + void set(const std::string& cache_id, + const GVector& vector); + void remove(const std::string& cache_id); + bool contains(const std::string& cache_id, + GVector* irfs = NULL) const; +}; + + +/***********************************************************************//** + * @brief GResponseVectorCache class extension + ***************************************************************************/ +%extend GResponseVectorCache { + GResponseVectorCache copy() { + return (*self); + } +}; diff --git a/pyext/gammalib/obs.i b/pyext/gammalib/obs.i index 50b3773..439f3d6 100644 --- a/pyext/gammalib/obs.i +++ b/pyext/gammalib/obs.i @@ -239,6 +239,7 @@ %include "GPhases.i" %include "GResponse.i" %include "GResponseCache.i" +%include "GResponseVectorCache.i" %include "GPhoton.i" %include "GPhotons.i" %include "GSource.i" diff --git a/src/obs/GResponse.cpp b/src/obs/GResponse.cpp index 80c587e..7c140e5 100644 --- a/src/obs/GResponse.cpp +++ b/src/obs/GResponse.cpp @@ -569,7 +569,7 @@ double GResponse::irf_spatial(const GEvent& event, double irf = 0.0; // Set IRF value attributes - std::string name = obs.id() + "::" + source.name(); + std::string name = obs.id() + ":" + source.name(); const GInstDir& dir = event.dir(); const GEnergy& ereco = event.energy(); const GEnergy& etrue = source.energy(); @@ -675,29 +675,49 @@ GVector GResponse::irf_spatial(const GModelSky& model, // Continue only if model has a spatial component if (model.spatial() != NULL) { - // Compute IRF for spatial model - switch (model.spatial()->code()) { - case GMODEL_SPATIAL_POINT_SOURCE: - irfs = irf_ptsrc(model, obs, gradients); - break; - case GMODEL_SPATIAL_RADIAL: - irfs = irf_radial(model, obs, gradients); - break; - case GMODEL_SPATIAL_ELLIPTICAL: - irfs = irf_elliptical(model, obs, gradients); - break; - case GMODEL_SPATIAL_DIFFUSE: - irfs = irf_diffuse(model, obs, gradients); - break; - case GMODEL_SPATIAL_COMPOSITE: - irfs = irf_composite(model, obs, gradients); - break; - default: - break; - } + // Set IRF cache identifier + std::string cache_id = obs.id() + ":" + model.name(); + + // Signal if spatial model has free parameters + bool has_free_pars = model.spatial()->has_free_pars(); + + // 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_vector_cache.contains(cache_id, &irfs)) { + + // Compute IRF for spatial model + switch (model.spatial()->code()) { + case GMODEL_SPATIAL_POINT_SOURCE: + irfs = irf_ptsrc(model, obs, gradients); + break; + case GMODEL_SPATIAL_RADIAL: + irfs = irf_radial(model, obs, gradients); + break; + case GMODEL_SPATIAL_ELLIPTICAL: + irfs = irf_elliptical(model, obs, gradients); + break; + case GMODEL_SPATIAL_DIFFUSE: + irfs = irf_diffuse(model, obs, gradients); + break; + case GMODEL_SPATIAL_COMPOSITE: + irfs = irf_composite(model, obs, gradients); + 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_vector_cache.set(cache_id, irfs); + } + + } // endif: computed spatial model + } // endif: model had spatial component // Return IRF values return irfs; @@ -745,6 +765,7 @@ void GResponse::init_members(void) // Clear cache m_irf_cache.clear(); m_nroi_cache.clear(); + m_irf_vector_cache.clear(); // Return return; @@ -768,6 +789,7 @@ void GResponse::copy_members(const GResponse& rsp) m_irf_diffuse_resolution = rsp.m_irf_diffuse_resolution; m_irf_cache = rsp.m_irf_cache; m_nroi_cache = rsp.m_nroi_cache; + m_irf_vector_cache = rsp.m_irf_vector_cache; // Return return; diff --git a/src/obs/GResponseVectorCache.cpp b/src/obs/GResponseVectorCache.cpp new file mode 100644 index 0000000..a190aa6 --- /dev/null +++ b/src/obs/GResponseVectorCache.cpp @@ -0,0 +1,480 @@ +/*************************************************************************** + * GResponseVectorCache.cpp - Response vector cache class * + * ----------------------------------------------------------------------- * + * copyright (C) 2020 by Juergen Knoedlseder * + * ----------------------------------------------------------------------- * + * * + * This program is free software: you can redistribute it and/or modify * + * it under the terms of the GNU General Public License as published by * + * the Free Software Foundation, either version 3 of the License, or * + * (at your option) any later version. * + * * + * This program is distributed in the hope that it will be useful, * + * but WITHOUT ANY WARRANTY; without even the implied warranty of * + * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the * + * GNU General Public License for more details. * + * * + * You should have received a copy of the GNU General Public License * + * along with this program. If not, see . * + * * + ***************************************************************************/ +/** + * @file GResponseVectorCache.cpp + * @brief Response vector cache class implementation + * @author Juergen Knoedlseder + */ + +/* __ Includes ___________________________________________________________ */ +#ifdef HAVE_CONFIG_H +#include +#endif +#include "GTools.hpp" +#include "GVector.hpp" +#include "GResponseVectorCache.hpp" + +/* __ Method name definitions ____________________________________________ */ + +/* __ Macros _____________________________________________________________ */ + +/* __ Coding definitions _________________________________________________ */ + +/* __ Debug definitions __________________________________________________ */ +//#define G_DEBUG_VECTOR_CACHE //!< Debug vector cache + + +/*========================================================================== + = = + = Constructors/destructors = + = = + ==========================================================================*/ + +/***********************************************************************//** + * @brief Void constructor + ***************************************************************************/ +GResponseVectorCache::GResponseVectorCache(void) +{ + // Initialise class members + init_members(); + + // Return + return; +} + + +/***********************************************************************//** + * @brief Copy constructor + * + * @param[in] cache Response vector cache. + ***************************************************************************/ +GResponseVectorCache::GResponseVectorCache(const GResponseVectorCache& cache) +{ + // Initialise class members + init_members(); + + // Copy members + copy_members(cache); + + // Return + return; +} + + +/***********************************************************************//** + * @brief Destructor + ***************************************************************************/ +GResponseVectorCache::~GResponseVectorCache(void) +{ + // Free members + free_members(); + + // Return + return; +} + + +/*========================================================================== + = = + = Operators = + = = + ==========================================================================*/ + +/***********************************************************************//** + * @brief Assignment operator + * + * @param[in] cache Response vector cache. + * @return Response vector cache. + ***************************************************************************/ +GResponseVectorCache& GResponseVectorCache::operator=(const GResponseVectorCache& cache) +{ + // Execute only if object is not identical + if (this != &cache) { + + // Free members + free_members(); + + // Initialise private members + init_members(); + + // Copy members + copy_members(cache); + + } // endif: object was not identical + + // Return this object + return *this; +} + + +/*========================================================================== + = = + = Public methods = + = = + ==========================================================================*/ + +/***********************************************************************//** + * @brief Clear response vector cache + ***************************************************************************/ +void GResponseVectorCache::clear(void) +{ + // Free members + free_members(); + + // Initialise private members + init_members(); + + // Return + return; +} + + +/***********************************************************************//** + * @brief Clone response cache + * + * @return Pointer to deep copy of response cache. + ***************************************************************************/ +GResponseVectorCache* GResponseVectorCache::clone(void) const +{ + return new GResponseVectorCache(*this); +} + + +/***********************************************************************//** + * @brief Set cache value + * + * @param[in] cache_id Cache identifier. + * @param[in] vector Cache vector. + * + * Set cache vector for a given @p cache_id. + ***************************************************************************/ +void GResponseVectorCache::set(const std::string& cache_id, + const GVector& vector) +{ + // Debug vector cache + #if defined(G_DEBUG_VECTOR_CACHE) + std::cout << "GResponseVectorCache::set("; + std::cout << cache_id << "," << vector.size() << "):"; + #endif + + // Initialise pointers to values and indices + double* values = NULL; + int* indices = NULL; + + // Find cache index + int index = find_cache(cache_id); + + // Determine size of vector cache + int entries = vector.non_zeros(); + + // Debug vector cache + #if defined(G_DEBUG_VECTOR_CACHE) + if (index != -1) { + std::cout << " entry found at index " << index; + std::cout << " with " << m_cache_entries[index] << " elements."; + std::cout << " Require " << entries << " elements."; + std::cout << std::endl; + } + else { + std::cout << " no entry found."; + std::cout << " Require " << entries << " elements."; + std::cout << std::endl; + } + #endif + + // If cache identifier does not yet exist then create new cache entry + // and allocate memory to hold cache values and indices + if (index == -1) { + if (entries > 0) { + values = new double[entries]; + indices = new int[entries]; + } + m_cache_ids.push_back(cache_id); + m_cache_entries.push_back(entries); + m_cache_values.push_back(values); + m_cache_indices.push_back(indices); + } + + // ... otherwise update cache existing entry by allocating new memory + // to hold cache values and indices in case that the number of entries + // have changed + else { + if (m_cache_entries[index] != entries) { + if (entries > 0) { + values = new double[entries]; + indices = new int[entries]; + } + if (m_cache_values[index] != NULL) delete [] m_cache_values[index]; + if (m_cache_indices[index] != NULL) delete [] m_cache_indices[index]; + m_cache_entries[index] = entries; + m_cache_values[index] = values; + m_cache_indices[index] = indices; + } + else { + values = m_cache_values[index]; + indices = m_cache_indices[index]; + } + } + + // Set cache vector + for (int i = 0; i < vector.size(); ++i) { + if (vector[i] != 0.0) { + *values++ = vector[i]; + *indices++ = i; + } + } + + // Return + return; +} + + +/***********************************************************************//** + * @brief Remove cache + * + * @param[in] cache_id Cache identifier. + * + * Remove cache for a given @p cache_id. + ***************************************************************************/ +void GResponseVectorCache::remove(const std::string& cache_id) +{ + // Find cache index + int index = find_cache(cache_id); + + // If cache index is valid then remove corresponding cache + if (index != -1) { + + // Free memory + if (m_cache_values[index] != NULL) delete [] m_cache_values[index]; + if (m_cache_indices[index] != NULL) delete [] m_cache_indices[index]; + + // Erase entry + m_cache_ids.erase(m_cache_ids.begin() + index); + m_cache_entries.erase(m_cache_entries.begin() + index); + m_cache_values.erase(m_cache_values.begin() + index); + m_cache_indices.erase(m_cache_indices.begin() + index); + + } // endif: cache index was valid + + // Return + return; +} + + +/***********************************************************************//** + * @brief Check if cache contains a value for specific parameters + * + * @param[in] cache_id Cache identifier. + * @param[in,out] vector Pointer to cached vector (only if found). + * @return True if cached value was found, otherwise false. + * + * Check if the cache contains a value for a given @p name, reconstructed + * energy @p ereco, and true energy @p etrue. + * + * If the @p value pointer argument is not NULL, the method will return the + * cached value through this argument in case that the value exists. + ***************************************************************************/ +bool GResponseVectorCache::contains(const std::string& cache_id, + GVector* vector) const +{ + // Debug vector cache + #if defined(G_DEBUG_VECTOR_CACHE) + std::cout << "GResponseVectorCache::contains(" << cache_id; + if (vector == NULL) { + std::cout << ", NULL):"; + } + else { + std::cout << "," << vector->size() << "):"; + } + #endif + + // Find cache index + int index = find_cache(cache_id); + + // Set contains flag + bool contains = (index != -1); + + // Debug vector cache + #if defined(G_DEBUG_VECTOR_CACHE) + if (contains) { + std::cout << " entry found at index " << index; + std::cout << " with " << m_cache_entries[index] << " elements."; + std::cout << std::endl; + } + else { + std::cout << " no entry found." << std::endl; + } + #endif + + // If cache contains the identifier and a pointer to a vector was + // provided then fill the vector + if (contains && vector != NULL) { + + // Get pointer to values and indices + double* values = m_cache_values[index]; + int* indices = m_cache_indices[index]; + + // Fill vector + for (int i = 0; i < m_cache_entries[index]; ++i, ++values, ++indices) { + (*vector)[*indices] = *values; + } + + } // endif: filled vector + + // Return containment flag + return contains; +} + + +/***********************************************************************//** + * @brief Print response cache + * + * @param[in] chatter Chattiness. + * @return String containing response cache information. + ***************************************************************************/ +std::string GResponseVectorCache::print(const GChatter& chatter) const +{ + // Initialise result string + std::string result; + + // Continue only if chatter is not silent + if (chatter != SILENT) { + + // Append header + result.append("=== GResponseVectorCache ==="); + + // Append total cache size + result.append("\n"+gammalib::parformat("Number of cached vectors")); + result.append(gammalib::str(size())); + + // Append cache information + for (int i = 0; i < size(); ++i) { + result.append("\n"+gammalib::parformat("Identifier")); + result.append(m_cache_ids[i]); + result.append(" ("+gammalib::str(m_cache_entries[i])+" values"); + } + + } // endif: chatter was not silent + + // Return result + return result; +} + + +/*========================================================================== + = = + = Private methods = + = = + ==========================================================================*/ + +/***********************************************************************//** + * @brief Initialise class members + ***************************************************************************/ +void GResponseVectorCache::init_members(void) +{ + // Initialise members + m_cache_ids.clear(); + m_cache_entries.clear(); + m_cache_values.clear(); + m_cache_indices.clear(); + + // Return + return; +} + + +/***********************************************************************//** + * @brief Copy class members + * + * @param[in] cache Response cache. + ***************************************************************************/ +void GResponseVectorCache::copy_members(const GResponseVectorCache& cache) +{ + // Copy members + m_cache_ids = cache.m_cache_ids; + m_cache_entries = cache.m_cache_entries; + + // Copy values and indices + for (int i = 0; i < size(); ++i) { + double* values = NULL; + int* indices = NULL; + if (m_cache_entries[i] > 0) { + values = new double[m_cache_entries[i]]; + indices = new int[m_cache_entries[i]]; + for (int k = 0; k < m_cache_entries[i]; ++k) { + values[k] = cache.m_cache_values[i][k]; + values[k] = cache.m_cache_indices[i][k]; + } + } + m_cache_values.push_back(values); + m_cache_indices.push_back(indices); + } + + // Return + return; +} + + +/***********************************************************************//** + * @brief Delete class members + ***************************************************************************/ +void GResponseVectorCache::free_members(void) +{ + // Loop over entries + for (int i = 0; i < size(); ++i) { + if (m_cache_values[i] != NULL) delete [] m_cache_values[i]; + if (m_cache_indices[i] != NULL) delete [] m_cache_indices[i]; + m_cache_values[i] = NULL; + m_cache_indices[i] = NULL; + } + + // Return + return; +} + + +/***********************************************************************//** + * @brief Find cache + * + * @param[in] cache_id Cache identifier. + * @return Cache index (-1 if no cache was found) + * + * Find cache for a given @p cache_id. If no cache was found the method + * returns -1. + ***************************************************************************/ +int GResponseVectorCache::find_cache(const std::string& cache_id) const +{ + // Initialise cache index + int index = -1; + + // Loop over cache + for (int i = 0; i < m_cache_ids.size(); ++i) { + if (m_cache_ids[i] == cache_id) { + index = i; + break; + } + } + + // Return index + return index; +} + + diff --git a/src/obs/Makefile.am b/src/obs/Makefile.am index 0d5e78f..cc227ec 100644 --- a/src/obs/Makefile.am +++ b/src/obs/Makefile.am @@ -47,6 +47,7 @@ sources = GEnergies.cpp \ GPhases.cpp \ GResponse.cpp \ GResponseCache.cpp \ + GResponseVectorCache.cpp \ GInstDir.cpp \ GPhotons.cpp \ GPhoton.cpp \ diff --git a/test/test_GObservation.cpp b/test/test_GObservation.cpp index ee0ca5a..4513d12 100644 --- a/test/test_GObservation.cpp +++ b/test/test_GObservation.cpp @@ -75,6 +75,8 @@ void TestGObservation::set(void) "Test GPhotons class"); append(static_cast(&TestGObservation::test_response_cache), "Test GResponseCache class"); + append(static_cast(&TestGObservation::test_response_vector_cache), + "Test GResponseVectorCache class"); append(static_cast(&TestGObservation::test_observations_optimizer), "Test GObservations::likelihood class"); @@ -1623,6 +1625,84 @@ void TestGObservation::test_response_cache(void) /***********************************************************************//** + * @brief Test GResponseCache + ***************************************************************************/ +void TestGObservation::test_response_vector_cache(void) +{ + // Allocate cache + GResponseVectorCache cache; + + // Set a few energies + GVector vector1(10); + vector1[1] = 1.0; + vector1[2] = 2.0; + GVector vector2(10); + vector2[5] = 5.0; + vector2[7] = 7.0; + GVector vector3(12); + vector3[1] = 1.0; + vector3[11] = 11.0; + + // Test empty cache + test_assert(cache.is_empty(), "Test that void cache is empty"); + test_value(cache.size(), 0, "Test void cache size"); + + // Add one vector and check + cache.set("Crab", vector1); + test_assert(!cache.is_empty(), "Test that filled cache is not empty"); + test_value(cache.size(), 1, "Test filled cache size"); + GVector vector4(10); + bool flag = cache.contains("Crab", &vector4); + test_assert(flag, "Test that cache contains Crab vector."); + test_assert(vector4 == vector1, "Test cache values Crab."); + + // Add another vector and check + cache.set("Vela", vector2); + test_assert(!cache.is_empty(), "Test that filled cache is not empty"); + test_value(cache.size(), 2, "Test filled cache size"); + GVector vector5(10); + flag = cache.contains("Vela", &vector5); + test_assert(flag, "Test that cache contains Vela vector."); + test_assert(vector5 == vector2, "Test cache values for Vela."); + + // Replace Crab values (same dimension) + cache.set("Crab", vector2); + test_assert(!cache.is_empty(), "Test that filled cache is not empty"); + test_value(cache.size(), 2, "Test filled cache size"); + GVector vector6(10); + flag = cache.contains("Crab", &vector6); + test_assert(flag, "Test that cache contains Crab vector."); + test_assert(vector6 == vector2, "Test cache values for replaced Crab."); + + // Replace Crab values (different dimension) + cache.set("Crab", vector3); + test_assert(!cache.is_empty(), "Test that filled cache is not empty"); + test_value(cache.size(), 2, "Test filled cache size"); + GVector vector7(12); + flag = cache.contains("Crab", &vector7); + test_assert(flag, "Test that cache contains Crab vector."); + test_assert(vector7 == vector3, "Test cache values for extended Crab."); + + // Remove Crab cache + cache.remove("Crab"); + test_value(cache.size(), 1, "Test filled cache size"); + flag = cache.contains("Crab", &vector7); + test_assert(!flag, "Test that cache does not contain Crab vector."); + flag = cache.contains("Vela", &vector5); + test_assert(flag, "Test that cache contains Vela vector."); + test_assert(vector5 == vector2, "Test cache values for Vela."); + + // Now clear the cache and test that it is empty + cache.clear(); + test_assert(cache.is_empty(), "Test that emptied cache is empty"); + test_value(cache.size(), 0, "Test emptied cache size"); + + // Return + return; +} + + +/***********************************************************************//** * @brief Test GObservations::likelihood ***************************************************************************/ void TestGObservation::test_observations_optimizer(void) @@ -1666,7 +1746,7 @@ void TestGObservation::test_observations_optimizer(void) // Free the event list delete events; - + } // endfor: added observations to container // Add the model to the observation @@ -1822,7 +1902,7 @@ void TestOpenMP::test_observations_optimizer(const int& mode) // Delete events pointer delete events; - + } // Add the model to the observation diff --git a/test/test_GObservation.hpp b/test/test_GObservation.hpp index ec969a1..0100c52 100644 --- a/test/test_GObservation.hpp +++ b/test/test_GObservation.hpp @@ -57,6 +57,7 @@ public: void test_energy(void); void test_energies(void); void test_response_cache(void); + void test_response_vector_cache(void); void test_observations_optimizer(void); }; @@ -68,7 +69,7 @@ public: * @brief Test suite for OpenMP testing ***************************************************************************/ class TestOpenMP : public GTestSuite { - + public: // Constructors and destructors TestOpenMP(void) : GTestSuite() {} -- libgit2 0.21.4