Luigi Tibaldo / gammalib

Commit d3768de7b6230a91048b26ce3bb52a6164d216e6

Authored by Jürgen Knödlseder
1 parent 97940c94

Implement sparse vector response cache (#3403)

Inline Side-by-side
Showing 15 changed files with 836 additions and 31 deletions
ChangeLog
  View file @d3768de
1   -2020-10-17
  1 +2020-10-20
2 2  
3 3 * Version 2.0.0 released
4 4 ========================
5 5  
  6 + Implement sparse vector response cache (#3403)
6 7 Implement GMatrixSparse::row() and GMatrixSparse::add_to_row() methods
7 8 Integrate analytical gradients for energy dispersion computation (#3333)
8 9 Add GFunctions abstract base class (#3328)
... ...
1 1 New Features and Important Changes in GammaLib 2.0.0
2 2  
3   -17 October 2020
  3 +20 October 2020
4 4  
5 5  
6 6 1. Introduction
... ... @@ -14,6 +14,7 @@ made since the last release of GammaLib.
14 14 The following classes have been added:
15 15 - GFunctions
16 16 - GIntegrals
  17 +- GResponseVectorCache
17 18  
18 19 The following classes have been removed:
19 20 - none
... ... @@ -87,6 +88,10 @@ of functions (#3328).
87 88  
88 89 10. Observation module
89 90 ----------------------
  91 +The GResponseVectorCache class was added to handle caching of the vector response.
  92 +Only non-zero IRF values are stored by the class. The class was integrated in the
  93 +GResponse class (#3403).
  94 +
90 95 A vectorised response computation was implemented in GResponse that computes the
91 96 IRF for all events in a given observations. This allows for a considerable reduction
92 97 of the operations, in particular for binned or stacked analysis, where the code can
... ...
README.md
  View file @d3768de
1 1 GammaLib information
2 2 ====================
3   -* Version: 2.0.0.dev (17 October 2020)
  3 +* Version: 2.0.0.dev (20 October 2020)
4 4  
5 5 [![Build Status](https://cta-jenkins.irap.omp.eu/buildStatus/icon?job=gammalib-integrate-os)](https://cta-jenkins.irap.omp.eu/job/gammalib-integrate-os/)
6 6  
... ...
doc/source/admin/release_history/2.0.rst
  View file @d3768de
1 1 .. _2.0:
2 2  
3   -GammaLib 2.0 release (17 October 2020)
  3 +GammaLib 2.0 release (20 October 2020)
4 4 ======================================
5 5  
6 6 Overview
... ... @@ -31,6 +31,8 @@ Improvements
31 31 New features
32 32 ------------
33 33  
  34 +* [`3403 <https://cta-redmine.irap.omp.eu/issues/3403>`_] -
  35 + Implement sparse vector response cache
34 36 * [`3328 <https://cta-redmine.irap.omp.eu/issues/3328>`_] -
35 37 Add ``GFunctions`` abstract base class and ``GIntegrals`` class
36 38  
... ...
include/GResponse.hpp
  View file @d3768de
... ... @@ -33,6 +33,7 @@
33 33 #include "GFunction.hpp"
34 34 #include "GTime.hpp"
35 35 #include "GResponseCache.hpp"
  36 +#include "GResponseVectorCache.hpp"
36 37 #include "GFunctions.hpp"
37 38  
38 39 /* __ Forward declarations _______________________________________________ */
... ... @@ -319,8 +320,9 @@ protected:
319 320 // Cache for irf(GEvent&, GSource&, GObservation&) and
320 321 // nroi(GModelSky&, GEnergy&, GTime&, GEnergy&, GTime&, GObservation&)
321 322 // computations
322   - mutable GResponseCache m_irf_cache;
323   - mutable GResponseCache m_nroi_cache;
  323 + mutable GResponseCache m_irf_cache;
  324 + mutable GResponseCache m_nroi_cache;
  325 + mutable GResponseVectorCache m_irf_vector_cache;
324 326 };
325 327  
326 328 #endif /* GRESPONSE_HPP */
... ...
include/GResponseVectorCache.hpp 0 → 100644
  View file @d3768de
  1 +/***************************************************************************
  2 + * GResponseVectorCache.hpp - Response vector cache class *
  3 + * ----------------------------------------------------------------------- *
  4 + * copyright (C) 2020 by Juergen Knoedlseder *
  5 + * ----------------------------------------------------------------------- *
  6 + * *
  7 + * This program is free software: you can redistribute it and/or modify *
  8 + * it under the terms of the GNU General Public License as published by *
  9 + * the Free Software Foundation, either version 3 of the License, or *
  10 + * (at your option) any later version. *
  11 + * *
  12 + * This program is distributed in the hope that it will be useful, *
  13 + * but WITHOUT ANY WARRANTY; without even the implied warranty of *
  14 + * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the *
  15 + * GNU General Public License for more details. *
  16 + * *
  17 + * You should have received a copy of the GNU General Public License *
  18 + * along with this program. If not, see <http://www.gnu.org/licenses/>. *
  19 + * *
  20 + ***************************************************************************/
  21 +/**
  22 + * @file GResponseVectorCache.hpp
  23 + * @brief Response vector cache class definition
  24 + * @author Juergen Knoedlseder
  25 + */
  26 +
  27 +#ifndef GRESPONSEVECTORCACHE_HPP
  28 +#define GRESPONSEVECTORCACHE_HPP
  29 +
  30 +/* __ Includes ___________________________________________________________ */
  31 +#include <vector>
  32 +#include <string>
  33 +#include "GBase.hpp"
  34 +
  35 +/* __ Forward declarations _______________________________________________ */
  36 +class GVector;
  37 +
  38 +/* __ Constants __________________________________________________________ */
  39 +
  40 +/* __ Typedefs ___________________________________________________________ */
  41 +
  42 +
  43 +/***********************************************************************//**
  44 + * @class GResponseVectorCache
  45 + *
  46 + * @brief Response vector cache class
  47 + *
  48 + * The class implements a vector cache for the Instrument Response Function
  49 + * values so that the values do not need to be recomputed each time but can
  50 + * be fetched from the cache.
  51 + ***************************************************************************/
  52 +class GResponseVectorCache : public GBase {
  53 +
  54 +public:
  55 + // Constructors and destructors
  56 + GResponseVectorCache(void);
  57 + GResponseVectorCache(const GResponseVectorCache& cache);
  58 + virtual ~GResponseVectorCache(void);
  59 +
  60 + // Operators
  61 + GResponseVectorCache& operator=(const GResponseVectorCache& cache);
  62 +
  63 + // Methods
  64 + void clear(void);
  65 + GResponseVectorCache* clone(void) const;
  66 + std::string classname(void) const;
  67 + bool is_empty(void) const;
  68 + int size(void) const;
  69 + void set(const std::string& cache_id,
  70 + const GVector& vector);
  71 + void remove(const std::string& cache_id);
  72 + bool contains(const std::string& cache_id,
  73 + GVector* irfs = NULL) const;
  74 + std::string print(const GChatter& chatter = NORMAL) const;
  75 +
  76 +protected:
  77 + // Protected methods
  78 + void init_members(void);
  79 + void copy_members(const GResponseVectorCache& cache);
  80 + void free_members(void);
  81 + int find_cache(const std::string& cache_id) const;
  82 +
  83 + // Protected members
  84 + std::vector<std::string> m_cache_ids;
  85 + std::vector<int> m_cache_entries;
  86 + std::vector<double*> m_cache_values;
  87 + std::vector<int*> m_cache_indices;
  88 +};
  89 +
  90 +
  91 +/***********************************************************************//**
  92 + * @brief Return class name
  93 + *
  94 + * @return String containing the class name ("GResponseVectorCache").
  95 + ***************************************************************************/
  96 +inline
  97 +std::string GResponseVectorCache::classname(void) const
  98 +{
  99 + return ("GResponseVectorCache");
  100 +}
  101 +
  102 +
  103 +/***********************************************************************//**
  104 + * @brief Checks whether the cache is empty
  105 + *
  106 + * @return True if cache is empty, false otherwise.
  107 + *
  108 + * Checks whether the response cache is empty.
  109 + ***************************************************************************/
  110 +inline
  111 +bool GResponseVectorCache::is_empty(void) const
  112 +{
  113 + return (m_cache_ids.empty());
  114 +}
  115 +
  116 +
  117 +/***********************************************************************//**
  118 + * @brief Remove cache for source
  119 + *
  120 + * @param[in] name Source name.
  121 + *
  122 + * Remove cache for source with @p name.
  123 + ***************************************************************************/
  124 +inline
  125 +int GResponseVectorCache::size(void) const
  126 +{
  127 + return (m_cache_ids.size());
  128 +}
  129 +
  130 +#endif /* GRESPONSEVECTORCACHE_HPP */
... ...
include/GammaLib.hpp
  View file @d3768de
... ... @@ -191,6 +191,7 @@
191 191 #include "GInstDir.hpp"
192 192 #include "GResponse.hpp"
193 193 #include "GResponseCache.hpp"
  194 +#include "GResponseVectorCache.hpp"
194 195 #include "GPhotons.hpp"
195 196 #include "GPhoton.hpp"
196 197 #include "GSource.hpp"
... ...
include/Makefile.am
  View file @d3768de
... ... @@ -151,6 +151,7 @@ pkginclude_HEADERS = GammaLib.hpp \
151 151 GInstDir.hpp \
152 152 GResponse.hpp \
153 153 GResponseCache.hpp \
  154 + GResponseVectorCache.hpp \
154 155 GPhotons.hpp \
155 156 GPhoton.hpp \
156 157 GSource.hpp \
... ...
pyext/GResponseVectorCache.i 0 → 100644
  View file @d3768de
  1 +/***************************************************************************
  2 + * GResponseVectorCache.i - Response vector cache class *
  3 + * ----------------------------------------------------------------------- *
  4 + * copyright (C) 2020 by Juergen Knoedlseder *
  5 + * ----------------------------------------------------------------------- *
  6 + * *
  7 + * This program is free software: you can redistribute it and/or modify *
  8 + * it under the terms of the GNU General Public License as published by *
  9 + * the Free Software Foundation, either version 3 of the License, or *
  10 + * (at your option) any later version. *
  11 + * *
  12 + * This program is distributed in the hope that it will be useful, *
  13 + * but WITHOUT ANY WARRANTY; without even the implied warranty of *
  14 + * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the *
  15 + * GNU General Public License for more details. *
  16 + * *
  17 + * You should have received a copy of the GNU General Public License *
  18 + * along with this program. If not, see <http://www.gnu.org/licenses/>. *
  19 + * *
  20 + ***************************************************************************/
  21 +/**
  22 + * @file GResponseVectorCache.i
  23 + * @brief Response vector cache class definition
  24 + * @author Juergen Knoedlseder
  25 + */
  26 +%{
  27 +/* Put headers and other declarations here that are needed for compilation */
  28 +#include "GResponseVectorCache.hpp"
  29 +%}
  30 +
  31 +/* __ Typemaps ___________________________________________________________ */
  32 +/*
  33 + * Typemap so that the contains() method returns two values, first the
  34 + * flag indicating whether the value exists or not, and second the
  35 + * value, i.e.
  36 + *
  37 + * flag, value = cache.contains('Crab', GEnergy(1.0,'TeV'), GEnergy(1.0,'TeV'))
  38 + *
  39 + */
  40 +%include "typemaps.i"
  41 +%apply GVector *OUTPUT { GVector* value };
  42 +
  43 +
  44 +/***********************************************************************//**
  45 + * @class GResponseVectorCache
  46 + *
  47 + * @brief Response cache class
  48 + ***************************************************************************/
  49 +class GResponseVectorCache : public GBase {
  50 +
  51 +public:
  52 + // Constructors and destructors
  53 + GResponseVectorCache(void);
  54 + GResponseVectorCache(const GResponseVectorCache& cache);
  55 + virtual ~GResponseVectorCache(void);
  56 +
  57 + // Methods
  58 + void clear(void);
  59 + GResponseVectorCache* clone(void) const;
  60 + std::string classname(void) const;
  61 + bool is_empty(void) const;
  62 + int size(void) const;
  63 + void set(const std::string& cache_id,
  64 + const GVector& vector);
  65 + void remove(const std::string& cache_id);
  66 + bool contains(const std::string& cache_id,
  67 + GVector* irfs = NULL) const;
  68 +};
  69 +
  70 +
  71 +/***********************************************************************//**
  72 + * @brief GResponseVectorCache class extension
  73 + ***************************************************************************/
  74 +%extend GResponseVectorCache {
  75 + GResponseVectorCache copy() {
  76 + return (*self);
  77 + }
  78 +};
... ...
pyext/gammalib/obs.i
  View file @d3768de
... ... @@ -239,6 +239,7 @@
239 239 %include "GPhases.i"
240 240 %include "GResponse.i"
241 241 %include "GResponseCache.i"
  242 +%include "GResponseVectorCache.i"
242 243 %include "GPhoton.i"
243 244 %include "GPhotons.i"
244 245 %include "GSource.i"
... ...
src/obs/GResponse.cpp
  View file @d3768de
... ... @@ -569,7 +569,7 @@ double GResponse::irf_spatial(const GEvent&amp; event,
569 569 double irf = 0.0;
570 570  
571 571 // Set IRF value attributes
572   - std::string name = obs.id() + "::" + source.name();
  572 + std::string name = obs.id() + ":" + source.name();
573 573 const GInstDir& dir = event.dir();
574 574 const GEnergy& ereco = event.energy();
575 575 const GEnergy& etrue = source.energy();
... ... @@ -675,29 +675,49 @@ GVector GResponse::irf_spatial(const GModelSky&amp; model,
675 675 // Continue only if model has a spatial component
676 676 if (model.spatial() != NULL) {
677 677  
678   - // Compute IRF for spatial model
679   - switch (model.spatial()->code()) {
680   - case GMODEL_SPATIAL_POINT_SOURCE:
681   - irfs = irf_ptsrc(model, obs, gradients);
682   - break;
683   - case GMODEL_SPATIAL_RADIAL:
684   - irfs = irf_radial(model, obs, gradients);
685   - break;
686   - case GMODEL_SPATIAL_ELLIPTICAL:
687   - irfs = irf_elliptical(model, obs, gradients);
688   - break;
689   - case GMODEL_SPATIAL_DIFFUSE:
690   - irfs = irf_diffuse(model, obs, gradients);
691   - break;
692   - case GMODEL_SPATIAL_COMPOSITE:
693   - irfs = irf_composite(model, obs, gradients);
694   - break;
695   - default:
696   - break;
697   - }
  678 + // Set IRF cache identifier
  679 + std::string cache_id = obs.id() + ":" + model.name();
  680 +
  681 + // Signal if spatial model has free parameters
  682 + bool has_free_pars = model.spatial()->has_free_pars();
  683 +
  684 + // If the spatial model component has free parameters, or the response
  685 + // cache should not be used, or the cache does not contain the requested
  686 + // IRF value then compute the IRF value for the spatial model.
  687 + if (has_free_pars ||
  688 + !m_use_irf_cache ||
  689 + !m_irf_vector_cache.contains(cache_id, &irfs)) {
  690 +
  691 + // Compute IRF for spatial model
  692 + switch (model.spatial()->code()) {
  693 + case GMODEL_SPATIAL_POINT_SOURCE:
  694 + irfs = irf_ptsrc(model, obs, gradients);
  695 + break;
  696 + case GMODEL_SPATIAL_RADIAL:
  697 + irfs = irf_radial(model, obs, gradients);
  698 + break;
  699 + case GMODEL_SPATIAL_ELLIPTICAL:
  700 + irfs = irf_elliptical(model, obs, gradients);
  701 + break;
  702 + case GMODEL_SPATIAL_DIFFUSE:
  703 + irfs = irf_diffuse(model, obs, gradients);
  704 + break;
  705 + case GMODEL_SPATIAL_COMPOSITE:
  706 + irfs = irf_composite(model, obs, gradients);
  707 + break;
  708 + default:
  709 + break;
  710 + }
698 711  
699   - } // endif: computed spatial model
  712 + // If the spatial model has no free parameters and the response cache
  713 + // should be used then put the IRF value in the response cache.
  714 + if (!has_free_pars && m_use_irf_cache) {
  715 + m_irf_vector_cache.set(cache_id, irfs);
  716 + }
  717 +
  718 + } // endif: computed spatial model
700 719  
  720 + } // endif: model had spatial component
701 721  
702 722 // Return IRF values
703 723 return irfs;
... ... @@ -745,6 +765,7 @@ void GResponse::init_members(void)
745 765 // Clear cache
746 766 m_irf_cache.clear();
747 767 m_nroi_cache.clear();
  768 + m_irf_vector_cache.clear();
748 769  
749 770 // Return
750 771 return;
... ... @@ -768,6 +789,7 @@ void GResponse::copy_members(const GResponse&amp; rsp)
768 789 m_irf_diffuse_resolution = rsp.m_irf_diffuse_resolution;
769 790 m_irf_cache = rsp.m_irf_cache;
770 791 m_nroi_cache = rsp.m_nroi_cache;
  792 + m_irf_vector_cache = rsp.m_irf_vector_cache;
771 793  
772 794 // Return
773 795 return;
... ...
src/obs/GResponseVectorCache.cpp 0 → 100644
  View file @d3768de
  1 +/***************************************************************************
  2 + * GResponseVectorCache.cpp - Response vector cache class *
  3 + * ----------------------------------------------------------------------- *
  4 + * copyright (C) 2020 by Juergen Knoedlseder *
  5 + * ----------------------------------------------------------------------- *
  6 + * *
  7 + * This program is free software: you can redistribute it and/or modify *
  8 + * it under the terms of the GNU General Public License as published by *
  9 + * the Free Software Foundation, either version 3 of the License, or *
  10 + * (at your option) any later version. *
  11 + * *
  12 + * This program is distributed in the hope that it will be useful, *
  13 + * but WITHOUT ANY WARRANTY; without even the implied warranty of *
  14 + * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the *
  15 + * GNU General Public License for more details. *
  16 + * *
  17 + * You should have received a copy of the GNU General Public License *
  18 + * along with this program. If not, see <http://www.gnu.org/licenses/>. *
  19 + * *
  20 + ***************************************************************************/
  21 +/**
  22 + * @file GResponseVectorCache.cpp
  23 + * @brief Response vector cache class implementation
  24 + * @author Juergen Knoedlseder
  25 + */
  26 +
  27 +/* __ Includes ___________________________________________________________ */
  28 +#ifdef HAVE_CONFIG_H
  29 +#include <config.h>
  30 +#endif
  31 +#include "GTools.hpp"
  32 +#include "GVector.hpp"
  33 +#include "GResponseVectorCache.hpp"
  34 +
  35 +/* __ Method name definitions ____________________________________________ */
  36 +
  37 +/* __ Macros _____________________________________________________________ */
  38 +
  39 +/* __ Coding definitions _________________________________________________ */
  40 +
  41 +/* __ Debug definitions __________________________________________________ */
  42 +//#define G_DEBUG_VECTOR_CACHE //!< Debug vector cache
  43 +
  44 +
  45 +/*==========================================================================
  46 + = =
  47 + = Constructors/destructors =
  48 + = =
  49 + ==========================================================================*/
  50 +
  51 +/***********************************************************************//**
  52 + * @brief Void constructor
  53 + ***************************************************************************/
  54 +GResponseVectorCache::GResponseVectorCache(void)
  55 +{
  56 + // Initialise class members
  57 + init_members();
  58 +
  59 + // Return
  60 + return;
  61 +}
  62 +
  63 +
  64 +/***********************************************************************//**
  65 + * @brief Copy constructor
  66 + *
  67 + * @param[in] cache Response vector cache.
  68 + ***************************************************************************/
  69 +GResponseVectorCache::GResponseVectorCache(const GResponseVectorCache& cache)
  70 +{
  71 + // Initialise class members
  72 + init_members();
  73 +
  74 + // Copy members
  75 + copy_members(cache);
  76 +
  77 + // Return
  78 + return;
  79 +}
  80 +
  81 +
  82 +/***********************************************************************//**
  83 + * @brief Destructor
  84 + ***************************************************************************/
  85 +GResponseVectorCache::~GResponseVectorCache(void)
  86 +{
  87 + // Free members
  88 + free_members();
  89 +
  90 + // Return
  91 + return;
  92 +}
  93 +
  94 +
  95 +/*==========================================================================
  96 + = =
  97 + = Operators =
  98 + = =
  99 + ==========================================================================*/
  100 +
  101 +/***********************************************************************//**
  102 + * @brief Assignment operator
  103 + *
  104 + * @param[in] cache Response vector cache.
  105 + * @return Response vector cache.
  106 + ***************************************************************************/
  107 +GResponseVectorCache& GResponseVectorCache::operator=(const GResponseVectorCache& cache)
  108 +{
  109 + // Execute only if object is not identical
  110 + if (this != &cache) {
  111 +
  112 + // Free members
  113 + free_members();
  114 +
  115 + // Initialise private members
  116 + init_members();
  117 +
  118 + // Copy members
  119 + copy_members(cache);
  120 +
  121 + } // endif: object was not identical
  122 +
  123 + // Return this object
  124 + return *this;
  125 +}
  126 +
  127 +
  128 +/*==========================================================================
  129 + = =
  130 + = Public methods =
  131 + = =
  132 + ==========================================================================*/
  133 +
  134 +/***********************************************************************//**
  135 + * @brief Clear response vector cache
  136 + ***************************************************************************/
  137 +void GResponseVectorCache::clear(void)
  138 +{
  139 + // Free members
  140 + free_members();
  141 +
  142 + // Initialise private members
  143 + init_members();
  144 +
  145 + // Return
  146 + return;
  147 +}
  148 +
  149 +
  150 +/***********************************************************************//**
  151 + * @brief Clone response cache
  152 + *
  153 + * @return Pointer to deep copy of response cache.
  154 + ***************************************************************************/
  155 +GResponseVectorCache* GResponseVectorCache::clone(void) const
  156 +{
  157 + return new GResponseVectorCache(*this);
  158 +}
  159 +
  160 +
  161 +/***********************************************************************//**
  162 + * @brief Set cache value
  163 + *
  164 + * @param[in] cache_id Cache identifier.
  165 + * @param[in] vector Cache vector.
  166 + *
  167 + * Set cache vector for a given @p cache_id.
  168 + ***************************************************************************/
  169 +void GResponseVectorCache::set(const std::string& cache_id,
  170 + const GVector& vector)
  171 +{
  172 + // Debug vector cache
  173 + #if defined(G_DEBUG_VECTOR_CACHE)
  174 + std::cout << "GResponseVectorCache::set(";
  175 + std::cout << cache_id << "," << vector.size() << "):";
  176 + #endif
  177 +
  178 + // Initialise pointers to values and indices
  179 + double* values = NULL;
  180 + int* indices = NULL;
  181 +
  182 + // Find cache index
  183 + int index = find_cache(cache_id);
  184 +
  185 + // Determine size of vector cache
  186 + int entries = vector.non_zeros();
  187 +
  188 + // Debug vector cache
  189 + #if defined(G_DEBUG_VECTOR_CACHE)
  190 + if (index != -1) {
  191 + std::cout << " entry found at index " << index;
  192 + std::cout << " with " << m_cache_entries[index] << " elements.";
  193 + std::cout << " Require " << entries << " elements.";
  194 + std::cout << std::endl;
  195 + }
  196 + else {
  197 + std::cout << " no entry found.";
  198 + std::cout << " Require " << entries << " elements.";
  199 + std::cout << std::endl;
  200 + }
  201 + #endif
  202 +
  203 + // If cache identifier does not yet exist then create new cache entry
  204 + // and allocate memory to hold cache values and indices
  205 + if (index == -1) {
  206 + if (entries > 0) {
  207 + values = new double[entries];
  208 + indices = new int[entries];
  209 + }
  210 + m_cache_ids.push_back(cache_id);
  211 + m_cache_entries.push_back(entries);
  212 + m_cache_values.push_back(values);
  213 + m_cache_indices.push_back(indices);
  214 + }
  215 +
  216 + // ... otherwise update cache existing entry by allocating new memory
  217 + // to hold cache values and indices in case that the number of entries
  218 + // have changed
  219 + else {
  220 + if (m_cache_entries[index] != entries) {
  221 + if (entries > 0) {
  222 + values = new double[entries];
  223 + indices = new int[entries];
  224 + }
  225 + if (m_cache_values[index] != NULL) delete [] m_cache_values[index];
  226 + if (m_cache_indices[index] != NULL) delete [] m_cache_indices[index];
  227 + m_cache_entries[index] = entries;
  228 + m_cache_values[index] = values;
  229 + m_cache_indices[index] = indices;
  230 + }
  231 + else {
  232 + values = m_cache_values[index];
  233 + indices = m_cache_indices[index];
  234 + }
  235 + }
  236 +
  237 + // Set cache vector
  238 + for (int i = 0; i < vector.size(); ++i) {
  239 + if (vector[i] != 0.0) {
  240 + *values++ = vector[i];
  241 + *indices++ = i;
  242 + }
  243 + }
  244 +
  245 + // Return
  246 + return;
  247 +}
  248 +
  249 +
  250 +/***********************************************************************//**
  251 + * @brief Remove cache
  252 + *
  253 + * @param[in] cache_id Cache identifier.
  254 + *
  255 + * Remove cache for a given @p cache_id.
  256 + ***************************************************************************/
  257 +void GResponseVectorCache::remove(const std::string& cache_id)
  258 +{
  259 + // Find cache index
  260 + int index = find_cache(cache_id);
  261 +
  262 + // If cache index is valid then remove corresponding cache
  263 + if (index != -1) {
  264 +
  265 + // Free memory
  266 + if (m_cache_values[index] != NULL) delete [] m_cache_values[index];
  267 + if (m_cache_indices[index] != NULL) delete [] m_cache_indices[index];
  268 +
  269 + // Erase entry
  270 + m_cache_ids.erase(m_cache_ids.begin() + index);
  271 + m_cache_entries.erase(m_cache_entries.begin() + index);
  272 + m_cache_values.erase(m_cache_values.begin() + index);
  273 + m_cache_indices.erase(m_cache_indices.begin() + index);
  274 +
  275 + } // endif: cache index was valid
  276 +
  277 + // Return
  278 + return;
  279 +}
  280 +
  281 +
  282 +/***********************************************************************//**
  283 + * @brief Check if cache contains a value for specific parameters
  284 + *
  285 + * @param[in] cache_id Cache identifier.
  286 + * @param[in,out] vector Pointer to cached vector (only if found).
  287 + * @return True if cached value was found, otherwise false.
  288 + *
  289 + * Check if the cache contains a value for a given @p name, reconstructed
  290 + * energy @p ereco, and true energy @p etrue.
  291 + *
  292 + * If the @p value pointer argument is not NULL, the method will return the
  293 + * cached value through this argument in case that the value exists.
  294 + ***************************************************************************/
  295 +bool GResponseVectorCache::contains(const std::string& cache_id,
  296 + GVector* vector) const
  297 +{
  298 + // Debug vector cache
  299 + #if defined(G_DEBUG_VECTOR_CACHE)
  300 + std::cout << "GResponseVectorCache::contains(" << cache_id;
  301 + if (vector == NULL) {
  302 + std::cout << ", NULL):";
  303 + }
  304 + else {
  305 + std::cout << "," << vector->size() << "):";
  306 + }
  307 + #endif
  308 +
  309 + // Find cache index
  310 + int index = find_cache(cache_id);
  311 +
  312 + // Set contains flag
  313 + bool contains = (index != -1);
  314 +
  315 + // Debug vector cache
  316 + #if defined(G_DEBUG_VECTOR_CACHE)
  317 + if (contains) {
  318 + std::cout << " entry found at index " << index;
  319 + std::cout << " with " << m_cache_entries[index] << " elements.";
  320 + std::cout << std::endl;
  321 + }
  322 + else {
  323 + std::cout << " no entry found." << std::endl;
  324 + }
  325 + #endif
  326 +
  327 + // If cache contains the identifier and a pointer to a vector was
  328 + // provided then fill the vector
  329 + if (contains && vector != NULL) {
  330 +
  331 + // Get pointer to values and indices
  332 + double* values = m_cache_values[index];
  333 + int* indices = m_cache_indices[index];
  334 +
  335 + // Fill vector
  336 + for (int i = 0; i < m_cache_entries[index]; ++i, ++values, ++indices) {
  337 + (*vector)[*indices] = *values;
  338 + }
  339 +
  340 + } // endif: filled vector
  341 +
  342 + // Return containment flag
  343 + return contains;
  344 +}
  345 +
  346 +
  347 +/***********************************************************************//**
  348 + * @brief Print response cache
  349 + *
  350 + * @param[in] chatter Chattiness.
  351 + * @return String containing response cache information.
  352 + ***************************************************************************/
  353 +std::string GResponseVectorCache::print(const GChatter& chatter) const
  354 +{
  355 + // Initialise result string
  356 + std::string result;
  357 +
  358 + // Continue only if chatter is not silent
  359 + if (chatter != SILENT) {
  360 +
  361 + // Append header
  362 + result.append("=== GResponseVectorCache ===");
  363 +
  364 + // Append total cache size
  365 + result.append("\n"+gammalib::parformat("Number of cached vectors"));
  366 + result.append(gammalib::str(size()));
  367 +
  368 + // Append cache information
  369 + for (int i = 0; i < size(); ++i) {
  370 + result.append("\n"+gammalib::parformat("Identifier"));
  371 + result.append(m_cache_ids[i]);
  372 + result.append(" ("+gammalib::str(m_cache_entries[i])+" values");
  373 + }
  374 +
  375 + } // endif: chatter was not silent
  376 +
  377 + // Return result
  378 + return result;
  379 +}
  380 +
  381 +
  382 +/*==========================================================================
  383 + = =
  384 + = Private methods =
  385 + = =
  386 + ==========================================================================*/
  387 +
  388 +/***********************************************************************//**
  389 + * @brief Initialise class members
  390 + ***************************************************************************/
  391 +void GResponseVectorCache::init_members(void)
  392 +{
  393 + // Initialise members
  394 + m_cache_ids.clear();
  395 + m_cache_entries.clear();
  396 + m_cache_values.clear();
  397 + m_cache_indices.clear();
  398 +
  399 + // Return
  400 + return;
  401 +}
  402 +
  403 +
  404 +/***********************************************************************//**
  405 + * @brief Copy class members
  406 + *
  407 + * @param[in] cache Response cache.
  408 + ***************************************************************************/
  409 +void GResponseVectorCache::copy_members(const GResponseVectorCache& cache)
  410 +{
  411 + // Copy members
  412 + m_cache_ids = cache.m_cache_ids;
  413 + m_cache_entries = cache.m_cache_entries;
  414 +
  415 + // Copy values and indices
  416 + for (int i = 0; i < size(); ++i) {
  417 + double* values = NULL;
  418 + int* indices = NULL;
  419 + if (m_cache_entries[i] > 0) {
  420 + values = new double[m_cache_entries[i]];
  421 + indices = new int[m_cache_entries[i]];
  422 + for (int k = 0; k < m_cache_entries[i]; ++k) {
  423 + values[k] = cache.m_cache_values[i][k];
  424 + values[k] = cache.m_cache_indices[i][k];
  425 + }
  426 + }
  427 + m_cache_values.push_back(values);
  428 + m_cache_indices.push_back(indices);
  429 + }
  430 +
  431 + // Return
  432 + return;
  433 +}
  434 +
  435 +
  436 +/***********************************************************************//**
  437 + * @brief Delete class members
  438 + ***************************************************************************/
  439 +void GResponseVectorCache::free_members(void)
  440 +{
  441 + // Loop over entries
  442 + for (int i = 0; i < size(); ++i) {
  443 + if (m_cache_values[i] != NULL) delete [] m_cache_values[i];
  444 + if (m_cache_indices[i] != NULL) delete [] m_cache_indices[i];
  445 + m_cache_values[i] = NULL;
  446 + m_cache_indices[i] = NULL;
  447 + }
  448 +
  449 + // Return
  450 + return;
  451 +}
  452 +
  453 +
  454 +/***********************************************************************//**
  455 + * @brief Find cache
  456 + *
  457 + * @param[in] cache_id Cache identifier.
  458 + * @return Cache index (-1 if no cache was found)
  459 + *
  460 + * Find cache for a given @p cache_id. If no cache was found the method
  461 + * returns -1.
  462 + ***************************************************************************/
  463 +int GResponseVectorCache::find_cache(const std::string& cache_id) const
  464 +{
  465 + // Initialise cache index
  466 + int index = -1;
  467 +
  468 + // Loop over cache
  469 + for (int i = 0; i < m_cache_ids.size(); ++i) {
  470 + if (m_cache_ids[i] == cache_id) {
  471 + index = i;
  472 + break;
  473 + }
  474 + }
  475 +
  476 + // Return index
  477 + return index;
  478 +}
  479 +
  480 +
... ...
src/obs/Makefile.am
  View file @d3768de
... ... @@ -47,6 +47,7 @@ sources = GEnergies.cpp \
47 47 GPhases.cpp \
48 48 GResponse.cpp \
49 49 GResponseCache.cpp \
  50 + GResponseVectorCache.cpp \
50 51 GInstDir.cpp \
51 52 GPhotons.cpp \
52 53 GPhoton.cpp \
... ...
test/test_GObservation.cpp
  View file @d3768de
... ... @@ -75,6 +75,8 @@ void TestGObservation::set(void)
75 75 "Test GPhotons class");
76 76 append(static_cast<pfunction>(&TestGObservation::test_response_cache),
77 77 "Test GResponseCache class");
  78 + append(static_cast<pfunction>(&TestGObservation::test_response_vector_cache),
  79 + "Test GResponseVectorCache class");
78 80 append(static_cast<pfunction>(&TestGObservation::test_observations_optimizer),
79 81 "Test GObservations::likelihood class");
80 82  
... ... @@ -1623,6 +1625,84 @@ void TestGObservation::test_response_cache(void)
1623 1625  
1624 1626  
1625 1627 /***********************************************************************//**
  1628 + * @brief Test GResponseCache
  1629 + ***************************************************************************/
  1630 +void TestGObservation::test_response_vector_cache(void)
  1631 +{
  1632 + // Allocate cache
  1633 + GResponseVectorCache cache;
  1634 +
  1635 + // Set a few energies
  1636 + GVector vector1(10);
  1637 + vector1[1] = 1.0;
  1638 + vector1[2] = 2.0;
  1639 + GVector vector2(10);
  1640 + vector2[5] = 5.0;
  1641 + vector2[7] = 7.0;
  1642 + GVector vector3(12);
  1643 + vector3[1] = 1.0;
  1644 + vector3[11] = 11.0;
  1645 +
  1646 + // Test empty cache
  1647 + test_assert(cache.is_empty(), "Test that void cache is empty");
  1648 + test_value(cache.size(), 0, "Test void cache size");
  1649 +
  1650 + // Add one vector and check
  1651 + cache.set("Crab", vector1);
  1652 + test_assert(!cache.is_empty(), "Test that filled cache is not empty");
  1653 + test_value(cache.size(), 1, "Test filled cache size");
  1654 + GVector vector4(10);
  1655 + bool flag = cache.contains("Crab", &vector4);
  1656 + test_assert(flag, "Test that cache contains Crab vector.");
  1657 + test_assert(vector4 == vector1, "Test cache values Crab.");
  1658 +
  1659 + // Add another vector and check
  1660 + cache.set("Vela", vector2);
  1661 + test_assert(!cache.is_empty(), "Test that filled cache is not empty");
  1662 + test_value(cache.size(), 2, "Test filled cache size");
  1663 + GVector vector5(10);
  1664 + flag = cache.contains("Vela", &vector5);
  1665 + test_assert(flag, "Test that cache contains Vela vector.");
  1666 + test_assert(vector5 == vector2, "Test cache values for Vela.");
  1667 +
  1668 + // Replace Crab values (same dimension)
  1669 + cache.set("Crab", vector2);
  1670 + test_assert(!cache.is_empty(), "Test that filled cache is not empty");
  1671 + test_value(cache.size(), 2, "Test filled cache size");
  1672 + GVector vector6(10);
  1673 + flag = cache.contains("Crab", &vector6);
  1674 + test_assert(flag, "Test that cache contains Crab vector.");
  1675 + test_assert(vector6 == vector2, "Test cache values for replaced Crab.");
  1676 +
  1677 + // Replace Crab values (different dimension)
  1678 + cache.set("Crab", vector3);
  1679 + test_assert(!cache.is_empty(), "Test that filled cache is not empty");
  1680 + test_value(cache.size(), 2, "Test filled cache size");
  1681 + GVector vector7(12);
  1682 + flag = cache.contains("Crab", &vector7);
  1683 + test_assert(flag, "Test that cache contains Crab vector.");
  1684 + test_assert(vector7 == vector3, "Test cache values for extended Crab.");
  1685 +
  1686 + // Remove Crab cache
  1687 + cache.remove("Crab");
  1688 + test_value(cache.size(), 1, "Test filled cache size");
  1689 + flag = cache.contains("Crab", &vector7);
  1690 + test_assert(!flag, "Test that cache does not contain Crab vector.");
  1691 + flag = cache.contains("Vela", &vector5);
  1692 + test_assert(flag, "Test that cache contains Vela vector.");
  1693 + test_assert(vector5 == vector2, "Test cache values for Vela.");
  1694 +
  1695 + // Now clear the cache and test that it is empty
  1696 + cache.clear();
  1697 + test_assert(cache.is_empty(), "Test that emptied cache is empty");
  1698 + test_value(cache.size(), 0, "Test emptied cache size");
  1699 +
  1700 + // Return
  1701 + return;
  1702 +}
  1703 +
  1704 +
  1705 +/***********************************************************************//**
1626 1706 * @brief Test GObservations::likelihood
1627 1707 ***************************************************************************/
1628 1708 void TestGObservation::test_observations_optimizer(void)
... ... @@ -1666,7 +1746,7 @@ void TestGObservation::test_observations_optimizer(void)
1666 1746  
1667 1747 // Free the event list
1668 1748 delete events;
1669   -
  1749 +
1670 1750 } // endfor: added observations to container
1671 1751  
1672 1752 // Add the model to the observation
... ... @@ -1822,7 +1902,7 @@ void TestOpenMP::test_observations_optimizer(const int&amp; mode)
1822 1902  
1823 1903 // Delete events pointer
1824 1904 delete events;
1825   -
  1905 +
1826 1906 }
1827 1907  
1828 1908 // Add the model to the observation
... ...
test/test_GObservation.hpp
  View file @d3768de
... ... @@ -57,6 +57,7 @@ public:
57 57 void test_energy(void);
58 58 void test_energies(void);
59 59 void test_response_cache(void);
  60 + void test_response_vector_cache(void);
60 61 void test_observations_optimizer(void);
61 62 };
62 63  
... ... @@ -68,7 +69,7 @@ public:
68 69 * @brief Test suite for OpenMP testing
69 70 ***************************************************************************/
70 71 class TestOpenMP : public GTestSuite {
71   -
  72 +
72 73 public:
73 74 // Constructors and destructors
74 75 TestOpenMP(void) : GTestSuite() {}
... ...