Commit 5d30b7bb747d87a6931fc21ce99fe606c9d6f3a2

Authored by Jürgen Knödlseder
1 parent 260282cf

Add GModelTemporalFunc class (#1914)

The GModelTemporalFunc implements a temporal file function for the simulation of transient sources.
ChangeLog
1   -2017-01-24
  1 +2017-01-25
2 2  
3 3 * Version 1.2.0 released
4 4 =========================
5 5  
  6 + Add GModelTemporalFunc call (#1914)
6 7 Add GSkyMap FITS HDU constructor
7 8 Add leap second for 1-1-2017 (#1887)
8 9 Add region() method to GModelSpatial (#1764)
... ...
1 1 New Features and Important Changes in GammaLib 1.2.0
2 2  
3   -24 January 2017
  3 +25 January 2017
4 4  
5 5  
6 6 1. Introduction
... ... @@ -17,6 +17,7 @@ The following classes have been added:
17 17 - GModelSpectralMultiplicative
18 18 - GModelSpectralComposite
19 19 - GModelSpatialComposite
  20 +- GModelTemporalFunc
20 21 - GVOTable
21 22 - GNdarray
22 23 - GFft
... ... @@ -169,6 +170,11 @@ derived classes that returns a pointer to a sky region that bounds the
169 170 spatial model. For the moment, only circular sky regions have been
170 171 implemented (#1764).
171 172  
  173 +The GModelTemporalFunc class has been added to implement a temporal file
  174 +function for the simulation of transient sources. A unit test has been
  175 +added and the mc() and eval() methods have been validated using a dedicated
  176 +test script (#1914).
  177 +
172 178  
173 179 9. Numerics module
174 180 ------------------
... ...
include/GModelTemporalFunc.hpp 0 → 100644
  1 +/***************************************************************************
  2 + * GModelTemporalFunc.hpp - Temporal file function model class *
  3 + * ----------------------------------------------------------------------- *
  4 + * copyright (C) 2017 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 GModelTemporalFunc.hpp
  23 + * @brief File function temporal model class interface definition
  24 + * @author Juergen Knoedlseder
  25 + */
  26 +
  27 +#ifndef GMODELTEMPORALFUNC_HPP
  28 +#define GMODELTEMPORALFUNC_HPP
  29 +
  30 +/* __ Includes ___________________________________________________________ */
  31 +#include <vector>
  32 +#include <string>
  33 +#include "GModelPar.hpp"
  34 +#include "GModelTemporal.hpp"
  35 +#include "GNodeArray.hpp"
  36 +#include "GFilename.hpp"
  37 +#include "GTimeReference.hpp"
  38 +#include "GTime.hpp"
  39 +
  40 +/* __ Forward declarations _______________________________________________ */
  41 +class GRan;
  42 +class GXmlElement;
  43 +
  44 +
  45 +/***********************************************************************//**
  46 + * @class GModelTemporalFunc
  47 + *
  48 + * @brief File function temporal model class
  49 + *
  50 + * This class implements a light curve defined by nodes given specified in
  51 + * a FITS file. The model is defined by
  52 + *
  53 + * \f[
  54 + * S_{\rm t}(t) = r(t) \times {\tt m\_norm}
  55 + * \f]
  56 + *
  57 + * where
  58 + * \f$r(t)$\f is the rate defined by linear interpolation between the nodes
  59 + * in a FITS file, and
  60 + * \f${\tt m\_norm}$\f is a normalisation constant.
  61 + ***************************************************************************/
  62 +class GModelTemporalFunc : public GModelTemporal {
  63 +
  64 +public:
  65 + // Constructors and destructors
  66 + GModelTemporalFunc(void);
  67 + explicit GModelTemporalFunc(const GXmlElement& xml);
  68 + GModelTemporalFunc(const GFilename& filename,
  69 + const double& norm = 1.0);
  70 + GModelTemporalFunc(const GModelTemporalFunc& model);
  71 + virtual ~GModelTemporalFunc(void);
  72 +
  73 + // Operators
  74 + virtual GModelTemporalFunc& operator=(const GModelTemporalFunc& model);
  75 +
  76 + // Implemented virtual base class methods
  77 + virtual void clear(void);
  78 + virtual GModelTemporalFunc* clone(void) const;
  79 + virtual std::string classname(void) const;
  80 + virtual std::string type(void) const;
  81 + virtual double eval(const GTime& srcTime,
  82 + const bool& gradients = false) const;
  83 + virtual GTimes mc(const double& rate, const GTime& tmin,
  84 + const GTime& tmax, GRan& ran) const;
  85 + virtual void read(const GXmlElement& xml);
  86 + virtual void write(GXmlElement& xml) const;
  87 + virtual std::string print(const GChatter& chatter = NORMAL) const;
  88 +
  89 + // Other methods
  90 + const GFilename& filename(void) const;
  91 + void filename(const GFilename& filename);
  92 + double norm(void) const;
  93 + void norm(const double& norm);
  94 +
  95 +protected:
  96 + // Protected methods
  97 + void init_members(void);
  98 + void copy_members(const GModelTemporalFunc& model);
  99 + void free_members(void);
  100 + void load_nodes(const GFilename& filename);
  101 + void mc_update(const GTime& tmin, const GTime& tmax) const;
  102 +
  103 + // Protected members
  104 + GModelPar m_norm; //!< Normalization factor
  105 + mutable GNodeArray m_nodes; //!< Time nodes of function
  106 + std::vector<double> m_values; //!< Function values at nodes
  107 + GFilename m_filename; //!< Name of file function
  108 + GTimeReference m_timeref; //!< Time reference
  109 + GTime m_tmin; //!< Minimum time of model
  110 + GTime m_tmax; //!< Maximum time of model
  111 +
  112 + // Cached members for MC
  113 + mutable GTime m_mc_tmin; //!< Minimum time
  114 + mutable GTime m_mc_tmax; //!< Maximum time
  115 + mutable double m_mc_rate; //!< Mean rate
  116 + mutable std::vector<double> m_mc_cum; //!< Cumulative distribution
  117 + mutable std::vector<double> m_mc_slope; //!< Slope of interval
  118 + mutable std::vector<double> m_mc_offset; //!< Offset of interval
  119 + mutable std::vector<double> m_mc_time; //!< Start time of interval
  120 + mutable std::vector<double> m_mc_dt; //!< Length of interval
  121 +};
  122 +
  123 +
  124 +/***********************************************************************//**
  125 + * @brief Return class name
  126 + *
  127 + * @return String containing the class name ("GModelTemporalFunc").
  128 + ***************************************************************************/
  129 +inline
  130 +std::string GModelTemporalFunc::classname(void) const
  131 +{
  132 + return ("GModelTemporalFunc");
  133 +}
  134 +
  135 +
  136 +/***********************************************************************//**
  137 + * @brief Return model type
  138 + *
  139 + * @return "FileFunction".
  140 + *
  141 + * Returns the type of the temporal model.
  142 + ***************************************************************************/
  143 +inline
  144 +std::string GModelTemporalFunc::type(void) const
  145 +{
  146 + return "FileFunction";
  147 +}
  148 +
  149 +
  150 +/***********************************************************************//**
  151 + * @brief Return normalization factor
  152 + *
  153 + * @return Normalization factor.
  154 + *
  155 + * Returns the normalization factor.
  156 + ***************************************************************************/
  157 +inline
  158 +double GModelTemporalFunc::norm(void) const
  159 +{
  160 + return (m_norm.value());
  161 +}
  162 +
  163 +
  164 +/***********************************************************************//**
  165 + * @brief Set normalization factor
  166 + *
  167 + * @param[in] norm Normalization factor.
  168 + *
  169 + * Sets the normalization factor.
  170 + ***************************************************************************/
  171 +inline
  172 +void GModelTemporalFunc::norm(const double& norm)
  173 +{
  174 + m_norm.value(norm);
  175 + return;
  176 +}
  177 +
  178 +
  179 +/***********************************************************************//**
  180 + * @brief Return file name
  181 + *
  182 + * @return Name of node file.
  183 + *
  184 + * Returns the name of the file function node file.
  185 + ***************************************************************************/
  186 +inline
  187 +const GFilename& GModelTemporalFunc::filename(void) const
  188 +{
  189 + return (m_filename);
  190 +}
  191 +
  192 +
  193 +/***********************************************************************//**
  194 + * @brief Loads nodes from node file and set filename
  195 + *
  196 + * @param[in] filename Node file name.
  197 + *
  198 + * Loads the nodes from a file function node file and sets the filename.
  199 + ***************************************************************************/
  200 +inline
  201 +void GModelTemporalFunc::filename(const GFilename& filename)
  202 +{
  203 + load_nodes(filename);
  204 + return;
  205 +}
  206 +
  207 +#endif /* GMODELTEMPORALFUNC_HPP */
... ...
include/GammaLib.hpp
1 1 /***************************************************************************
2 2 * GammaLib.hpp - Gamma-Ray Astronomy Library Header file *
3 3 * ----------------------------------------------------------------------- *
4   - * copyright (C) 2006-2016 by Juergen Knoedlseder *
  4 + * copyright (C) 2006-2017 by Juergen Knoedlseder *
5 5 * ----------------------------------------------------------------------- *
6 6 * *
7 7 * This program is free software: you can redistribute it and/or modify *
... ... @@ -228,5 +228,6 @@
228 228 #include "GModelTemporal.hpp"
229 229 #include "GModelTemporalRegistry.hpp"
230 230 #include "GModelTemporalConst.hpp"
  231 +#include "GModelTemporalFunc.hpp"
231 232  
232 233 #endif /* GAMMALIB_HPP */
... ...
include/Makefile.am
1 1 #############################################################################
2 2 # Makefile.am for include module #
3 3 # ------------------------------------------------------------------------- #
4   -# Copyright (C) 2016 Juergen Knoedlseder #
  4 +# Copyright (C) 2016-2017 Juergen Knoedlseder #
5 5 # ------------------------------------------------------------------------- #
6 6 # #
7 7 # This program is free software: you can redistribute it and/or modify #
... ... @@ -185,4 +185,5 @@ pkginclude_HEADERS = GammaLib.hpp \
185 185 GModelSpectralComposite.hpp \
186 186 GModelTemporal.hpp \
187 187 GModelTemporalRegistry.hpp \
188   - GModelTemporalConst.hpp
  188 + GModelTemporalConst.hpp \
  189 + GModelTemporalFunc.hpp
... ...
pyext/GModelTemporalFunc.i 0 → 100644
  1 +/***************************************************************************
  2 + * GModelTemporalFunc.i - Temporal file function model class *
  3 + * ----------------------------------------------------------------------- *
  4 + * copyright (C) 2017 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 GModelTemporalFunc.i
  23 + * @brief File function temporal model class interface definition
  24 + * @author Juergen Knoedlseder
  25 + */
  26 +%{
  27 +/* Put headers and other declarations here that are needed for compilation */
  28 +#include "GModelTemporalFunc.hpp"
  29 +%}
  30 +
  31 +
  32 +/***********************************************************************//**
  33 + * @class GModelTemporalFunc
  34 + *
  35 + * @brief File function temporal model class
  36 + ***************************************************************************/
  37 +class GModelTemporalFunc : public GModelTemporal {
  38 +
  39 +public:
  40 + // Constructors and destructors
  41 + GModelTemporalFunc(void);
  42 + explicit GModelTemporalFunc(const GXmlElement& xml);
  43 + GModelTemporalFunc(const GFilename& filename,
  44 + const double& norm = 1.0);
  45 + GModelTemporalFunc(const GModelTemporalFunc& model);
  46 + virtual ~GModelTemporalFunc(void);
  47 +
  48 + // Implemented virtual base class methods
  49 + virtual void clear(void);
  50 + virtual GModelTemporalFunc* clone(void) const;
  51 + virtual std::string classname(void) const;
  52 + virtual std::string type(void) const;
  53 + virtual double eval(const GTime& srcTime,
  54 + const bool& gradients = false) const;
  55 + virtual GTimes mc(const double& rate, const GTime& tmin,
  56 + const GTime& tmax, GRan& ran) const;
  57 + virtual void read(const GXmlElement& xml);
  58 + virtual void write(GXmlElement& xml) const;
  59 +
  60 + // Other methods
  61 + const GFilename& filename(void) const;
  62 + void filename(const GFilename& filename);
  63 + double norm(void) const;
  64 + void norm(const double& norm);
  65 +};
  66 +
  67 +
  68 +/***********************************************************************//**
  69 + * @brief GModelTemporalFunc class extension
  70 + ***************************************************************************/
  71 +%extend GModelTemporalFunc {
  72 + GModelTemporalFunc copy() {
  73 + return (*self);
  74 + }
  75 +};
... ...
pyext/gammalib/model.i
1 1 /***************************************************************************
2 2 * model.i - Model module *
3 3 * ----------------------------------------------------------------------- *
4   - * copyright (C) 2011-2016 by Juergen Knoedlseder *
  4 + * copyright (C) 2011-2017 by Juergen Knoedlseder *
5 5 * ----------------------------------------------------------------------- *
6 6 * *
7 7 * This program is free software: you can redistribute it and/or modify *
... ... @@ -160,3 +160,4 @@
160 160 %include "GModelTemporal.i"
161 161 %include "GModelTemporalRegistry.i"
162 162 %include "GModelTemporalConst.i"
  163 +%include "GModelTemporalFunc.i"
... ...
src/model/GModelTemporalFunc.cpp 0 → 100644
  1 +/***************************************************************************
  2 + * GModelTemporalFunc.cpp - Temporal file function model class *
  3 + * ----------------------------------------------------------------------- *
  4 + * copyright (C) 2017 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 GModelTemporalFunc.cpp
  23 + * @brief File function temporal model class interface 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 "GException.hpp"
  33 +#include "GModelTemporalFunc.hpp"
  34 +#include "GModelTemporalRegistry.hpp"
  35 +#include "GFits.hpp"
  36 +
  37 +/* __ Constants __________________________________________________________ */
  38 +
  39 +/* __ Globals ____________________________________________________________ */
  40 +const GModelTemporalFunc g_temporal_func_seed;
  41 +const GModelTemporalRegistry g_temporal_func_registry(&g_temporal_func_seed);
  42 +
  43 +/* __ Method name definitions ____________________________________________ */
  44 +#define G_READ "GModelTemporalFunc::read(GXmlElement&)"
  45 +#define G_WRITE "GModelTemporalFunc::write(GXmlElement&)"
  46 +#define G_LOAD_NODES "GModelTemporalFunc::load_nodes(GFilename&)"
  47 +
  48 +/* __ Macros _____________________________________________________________ */
  49 +
  50 +/* __ Coding definitions _________________________________________________ */
  51 +
  52 +/* __ Debug definitions __________________________________________________ */
  53 +
  54 +
  55 +/*==========================================================================
  56 + = =
  57 + = Constructors/destructors =
  58 + = =
  59 + ==========================================================================*/
  60 +
  61 +/***********************************************************************//**
  62 + * @brief Void constructor
  63 + ***************************************************************************/
  64 +GModelTemporalFunc::GModelTemporalFunc(void) : GModelTemporal()
  65 +{
  66 + // Initialise members
  67 + init_members();
  68 +
  69 + // Return
  70 + return;
  71 +}
  72 +
  73 +
  74 +/***********************************************************************//**
  75 + * @brief File constructor
  76 + *
  77 + * @param[in] filename File name of nodes.
  78 + * @param[in] norm Normalization factor.
  79 + *
  80 + * Constructs spectral file function model from a list of nodes that is found
  81 + * in the specified FITS file. See the load_nodes() method for more
  82 + * information about the expected structure of the file.
  83 + ***************************************************************************/
  84 +GModelTemporalFunc::GModelTemporalFunc(const GFilename& filename,
  85 + const double& norm) :
  86 + GModelTemporal()
  87 +{
  88 + // Initialise members
  89 + init_members();
  90 +
  91 + // Load nodes
  92 + load_nodes(filename);
  93 +
  94 + // Set normalization
  95 + m_norm.value(norm);
  96 +
  97 + // Return
  98 + return;
  99 +}
  100 +
  101 +
  102 +/***********************************************************************//**
  103 + * @brief XML constructor
  104 + *
  105 + * @param[in] xml XML element.
  106 + *
  107 + * Constructs constant temporal model by extracting information from an XML
  108 + * element. See the read() method for more information about the expected
  109 + * structure of the XML element.
  110 + ***************************************************************************/
  111 +GModelTemporalFunc::GModelTemporalFunc(const GXmlElement& xml) :
  112 + GModelTemporal()
  113 +{
  114 + // Initialise members
  115 + init_members();
  116 +
  117 + // Read information from XML element
  118 + read(xml);
  119 +
  120 + // Return
  121 + return;
  122 +}
  123 +
  124 +
  125 +/***********************************************************************//**
  126 + * @brief Copy constructor
  127 + *
  128 + * @param[in] model File function temporal model.
  129 + ***************************************************************************/
  130 +GModelTemporalFunc::GModelTemporalFunc(const GModelTemporalFunc& model) :
  131 + GModelTemporal(model)
  132 +{
  133 + // Initialise members
  134 + init_members();
  135 +
  136 + // Copy members
  137 + copy_members(model);
  138 +
  139 + // Return
  140 + return;
  141 +}
  142 +
  143 +
  144 +/***********************************************************************//**
  145 + * @brief Destructor
  146 + ***************************************************************************/
  147 +GModelTemporalFunc::~GModelTemporalFunc(void)
  148 +{
  149 + // Free members
  150 + free_members();
  151 +
  152 + // Return
  153 + return;
  154 +}
  155 +
  156 +
  157 +/*==========================================================================
  158 + = =
  159 + = Operators =
  160 + = =
  161 + ==========================================================================*/
  162 +
  163 +/***********************************************************************//**
  164 + * @brief Assignment operator
  165 + *
  166 + * @param[in] model File function temporal model.
  167 + * @return File function temporal model.
  168 + ***************************************************************************/
  169 +GModelTemporalFunc& GModelTemporalFunc::operator=(const GModelTemporalFunc& model)
  170 +{
  171 + // Execute only if object is not identical
  172 + if (this != &model) {
  173 +
  174 + // Copy base class members
  175 + this->GModelTemporal::operator=(model);
  176 +
  177 + // Free members
  178 + free_members();
  179 +
  180 + // Initialise members
  181 + init_members();
  182 +
  183 + // Copy members
  184 + copy_members(model);
  185 +
  186 + } // endif: object was not identical
  187 +
  188 + // Return
  189 + return *this;
  190 +}
  191 +
  192 +
  193 +/*==========================================================================
  194 + = =
  195 + = Public methods =
  196 + = =
  197 + ==========================================================================*/
  198 +
  199 +/***********************************************************************//**
  200 + * @brief Clear file function temporal model
  201 + ***************************************************************************/
  202 +void GModelTemporalFunc::clear(void)
  203 +{
  204 + // Free class members (base and derived classes, derived class first)
  205 + free_members();
  206 + this->GModelTemporal::free_members();
  207 +
  208 + // Initialise members
  209 + this->GModelTemporal::init_members();
  210 + init_members();
  211 +
  212 + // Return
  213 + return;
  214 +}
  215 +
  216 +
  217 +/***********************************************************************//**
  218 + * @brief Clone file function temporal model
  219 + *
  220 + * @return Pointer to deep copy of file function temporal model.
  221 + ***************************************************************************/
  222 +GModelTemporalFunc* GModelTemporalFunc::clone(void) const
  223 +{
  224 + // Clone constant temporal model
  225 + return new GModelTemporalFunc(*this);
  226 +}
  227 +
  228 +
  229 +/***********************************************************************//**
  230 + * @brief Evaluate function
  231 + *
  232 + * @param[in] srcTime True photon arrival time (not used).
  233 + * @param[in] gradients Compute gradients?
  234 + * @return Value of temporal model.
  235 + *
  236 + * Computes
  237 + *
  238 + * \f[
  239 + * S_{\rm t}(t) = r(t) \times {\tt m\_norm}
  240 + * \f]
  241 + *
  242 + * where
  243 + * \f$r(t)\f$ is the linear interpolation between the nodes of the file
  244 + * function and
  245 + * \f${\tt m\_norm}\f$ is the normalization constant.
  246 + *
  247 + * If the @p gradients flag is true the method will also evaluate the partial
  248 + * derivatives of the model with respect to the normalization parameter using
  249 + *
  250 + * \f[
  251 + * \frac{\delta S_{\rm t}(t)}{\delta {\tt m\_norm}} = r(t)
  252 + * \f]
  253 + ***************************************************************************/
  254 +double GModelTemporalFunc::eval(const GTime& srcTime,
  255 + const bool& gradients) const
  256 +{
  257 + // Initialise value
  258 + double value = 0.0;
  259 +
  260 + // Optionally initialise gradient
  261 + if (gradients) {
  262 + m_norm.factor_gradient(0.0);
  263 + }
  264 +
  265 + // Set value if the time is within the validity range
  266 + if (srcTime >= m_tmin && srcTime <= m_tmax) {
  267 +
  268 + // Interpolate file function
  269 + double func = m_nodes.interpolate(srcTime.convert(m_timeref), m_values);
  270 +
  271 + // Compute function value
  272 + value = m_norm.value() * func;
  273 +
  274 + // Optionally compute gradients
  275 + if (gradients) {
  276 +
  277 + // Compute partial derivatives of the parameter values
  278 + double g_norm = (m_norm.is_free()) ? m_norm.scale() * func : 0.0;
  279 +
  280 + // Set gradients
  281 + m_norm.factor_gradient(g_norm);
  282 +
  283 + } // endif: gradient computation was requested
  284 +
  285 + } // endif: time was within the validity range
  286 +
  287 + // Return
  288 + return value;
  289 +}
  290 +
  291 +
  292 +/***********************************************************************//**
  293 + * @brief Returns vector of random event times
  294 + *
  295 + * @param[in] rate Mean event rate (events per second).
  296 + * @param[in] tmin Minimum event time.
  297 + * @param[in] tmax Maximum event time.
  298 + * @param[in,out] ran Random number generator.
  299 + *
  300 + * This method returns a vector of random event times between @p tmin and
  301 + * @p tmax assuming a light curve specified by a file function.
  302 + ***************************************************************************/
  303 +GTimes GModelTemporalFunc::mc(const double& rate, const GTime& tmin,
  304 + const GTime& tmax, GRan& ran) const
  305 +{
  306 + // Allocates empty vector of times
  307 + GTimes times;
  308 +
  309 + // Update Monte Carlo cache
  310 + mc_update(tmin, tmax);
  311 +
  312 + // Continue only if mean rate is positive and cache is not empty
  313 + if (m_mc_rate > 0.0 && m_mc_cum.size() > 0) {
  314 +
  315 + // Compute mean event rate (in events per seconds)
  316 + double lambda = rate * norm() * m_mc_rate;
  317 +
  318 + // Compute number of times to be sampled
  319 + int ntimes = int(ran.poisson(lambda)+0.5);
  320 +
  321 + // Loop over number of times
  322 + for (int i = 0; i < ntimes; ++i) {
  323 +
  324 + // Determine in which bin we reside
  325 + int inx = 0;
  326 + if (m_mc_cum.size() > 1) {
  327 + double u = ran.uniform();
  328 + for (inx = m_mc_cum.size()-1; inx > 0; --inx) {
  329 + if (m_mc_cum[inx-1] <= u) {
  330 + break;
  331 + }
  332 + }
  333 + }
  334 +
  335 + // Get random time
  336 + double seconds;
  337 + if (m_mc_slope[inx] == 0.0) {
  338 + seconds = m_mc_dt[inx] * ran.uniform() + m_mc_time[inx];
  339 + }
  340 + else {
  341 + seconds = (std::sqrt(m_mc_offset[inx]*m_mc_offset[inx] +
  342 + 2.0 * m_mc_slope[inx] * ran.uniform()) -
  343 + m_mc_offset[inx]) / m_mc_slope[inx] + m_mc_time[inx];
  344 + }
  345 +
  346 + // Append random time
  347 + GTime time(seconds, m_timeref);
  348 + times.append(time);
  349 +
  350 + } // endfor: looped over times
  351 +
  352 + } // endif: cache was valid
  353 +
  354 + // Return vector of times
  355 + return times;
  356 +}
  357 +
  358 +
  359 +/***********************************************************************//**
  360 + * @brief Read model from XML element
  361 + *
  362 + * @param[in] xml XML element.
  363 + *
  364 + * Writes the temporal information from an XML element having the format
  365 + *
  366 + * <temporalModel type="FileFunction" file="..">
  367 + * <parameter name="Normalization" scale="1" value="1" min="0.1" max="10" free="1"/>
  368 + * </temporalModel>
  369 + ***************************************************************************/
  370 +void GModelTemporalFunc::read(const GXmlElement& xml)
  371 +{
  372 + // Get parameter pointers
  373 + const GXmlElement* norm = gammalib::xml_get_par(G_READ, xml, m_norm.name());
  374 +
  375 + // Read parameters
  376 + m_norm.read(*norm);
  377 +
  378 + // Load nodes from file
  379 + load_nodes(gammalib::xml_file_expand(xml, xml.attribute("file")));
  380 +
  381 + // Return
  382 + return;
  383 +}
  384 +
  385 +
  386 +/***********************************************************************//**
  387 + * @brief Write model into XML element
  388 + *
  389 + * @param[in] xml XML element.
  390 + *
  391 + * @exception GExpection::invalid_value
  392 + * Invalid XML format encountered.
  393 + *
  394 + * Writes the temporal information into an XML element in the format
  395 + *
  396 + * <temporalModel type="Constant">
  397 + * <parameter name="Normalization" scale="1" value="1" min="0.1" max="10" free="1"/>
  398 + * </temporalModel>
  399 + ***************************************************************************/
  400 +void GModelTemporalFunc::write(GXmlElement& xml) const
  401 +{
  402 + // Set model type
  403 + if (xml.attribute("type") == "") {
  404 + xml.attribute("type", type());
  405 + }
  406 +
  407 + // Verify model type
  408 + if (xml.attribute("type") != type()) {
  409 + std::string msg = "Temporal model of type "+xml.attribute("type")+
  410 + " encountered while \""+type()+"\" was expected.";
  411 + throw GException::invalid_value(G_WRITE, msg);
  412 + }
  413 +
  414 + // Get XML parameters
  415 + GXmlElement* norm = gammalib::xml_need_par(G_WRITE, xml, m_norm.name());
  416 +
  417 + // Write parameters
  418 + m_norm.write(*norm);
  419 +
  420 + // Set file attribute
  421 + xml.attribute("file", gammalib::xml_file_reduce(xml, m_filename));
  422 +
  423 + // Return
  424 + return;
  425 +}
  426 +
  427 +
  428 +/***********************************************************************//**
  429 + * @brief Print file function information
  430 + *
  431 + * @param[in] chatter Chattiness.
  432 + * @return String containing model information.
  433 + ***************************************************************************/
  434 +std::string GModelTemporalFunc::print(const GChatter& chatter) const
  435 +{
  436 + // Initialise result string
  437 + std::string result;
  438 +
  439 + // Continue only if chatter is not silent
  440 + if (chatter != SILENT) {
  441 +
  442 + // Append header
  443 + result.append("=== GModelTemporalFunc ===");
  444 +
  445 + // Append information
  446 + result.append("\n"+gammalib::parformat("Function file"));
  447 + result.append(m_filename.url());
  448 + result.append("\n"+gammalib::parformat("Number of parameters"));
  449 + result.append(gammalib::str(size()));
  450 + for (int i = 0; i < size(); ++i) {
  451 + result.append("\n"+m_pars[i]->print(chatter));
  452 + }
  453 +
  454 + // Append reference MJD and time range
  455 + result.append("\n"+gammalib::parformat("Reference MDJ"));
  456 + result.append(gammalib::str(m_timeref.mjdref()));
  457 + result.append("\n"+gammalib::parformat("Time range"));
  458 + result.append(gammalib::str(m_tmin.convert(m_timeref)));
  459 + result.append(" - ");
  460 + result.append(gammalib::str(m_tmax.convert(m_timeref)));
  461 + result.append(" "+m_timeref.timeunit());
  462 + result.append(" ("+m_timeref.timesys()+")");
  463 +
  464 + // Append node information
  465 + result.append("\n"+gammalib::parformat("Number of nodes"));
  466 + result.append(gammalib::str(m_nodes.size()));
  467 + for (int i = 0; i < m_nodes.size(); ++i) {
  468 + result.append("\n");
  469 + result.append(gammalib::parformat(" Time "+gammalib::str(m_nodes[i])+" s"));
  470 + result.append(gammalib::str(m_values[i]));
  471 + }
  472 +
  473 + } // endif: chatter was not silent
  474 +
  475 + // Return result
  476 + return result;
  477 +}
  478 +
  479 +
  480 +/*==========================================================================
  481 + = =
  482 + = Private methods =
  483 + = =
  484 + ==========================================================================*/
  485 +
  486 +/***********************************************************************//**
  487 + * @brief Initialise class members
  488 + ***************************************************************************/
  489 +void GModelTemporalFunc::init_members(void)
  490 +{
  491 + // Initialise normalisation parameter
  492 + m_norm.clear();
  493 + m_norm.name("Normalization");
  494 + m_norm.unit("(relative value)");
  495 + m_norm.scale(1.0);
  496 + m_norm.value(1.0);
  497 + m_norm.range(0.0,1000.0);
  498 + m_norm.fix();
  499 + m_norm.gradient(0.0);
  500 + m_norm.has_grad(true);
  501 +
  502 + // Set parameter pointer(s)
  503 + m_pars.clear();
  504 + m_pars.push_back(&m_norm);
  505 +
  506 + // Initialise members
  507 + m_nodes.clear();
  508 + m_values.clear();
  509 + m_filename.clear();
  510 + m_timeref.clear();
  511 + m_tmin.clear();
  512 + m_tmax.clear();
  513 +
  514 + // Initialise cache
  515 + m_mc_tmin.clear();
  516 + m_mc_tmax.clear();
  517 + m_mc_rate = 0.0;
  518 + m_mc_cum.clear();
  519 + m_mc_slope.clear();
  520 + m_mc_offset.clear();
  521 + m_mc_time.clear();
  522 + m_mc_dt.clear();
  523 +
  524 + // Return
  525 + return;
  526 +}
  527 +
  528 +
  529 +/***********************************************************************//**
  530 + * @brief Copy class members
  531 + *
  532 + * @param[in] model Constant temporal model
  533 + ***************************************************************************/
  534 +void GModelTemporalFunc::copy_members(const GModelTemporalFunc& model)
  535 +{
  536 + // Copy members
  537 + m_norm = model.m_norm;
  538 + m_nodes = model.m_nodes;
  539 + m_values = model.m_values;
  540 + m_filename = model.m_filename;
  541 + m_timeref = model.m_timeref;
  542 + m_tmin = model.m_tmin;
  543 + m_tmax = model.m_tmax;
  544 +
  545 + // Copy cache
  546 + m_mc_tmin = model.m_mc_tmin;
  547 + m_mc_tmax = model.m_mc_tmax;
  548 + m_mc_rate = model.m_mc_rate;
  549 + m_mc_cum = model.m_mc_cum;
  550 + m_mc_slope = model.m_mc_slope;
  551 + m_mc_offset = model.m_mc_offset;
  552 + m_mc_time = model.m_mc_time;
  553 + m_mc_dt = model.m_mc_dt;
  554 +
  555 + // Set parameter pointer(s)
  556 + m_pars.clear();
  557 + m_pars.push_back(&m_norm);
  558 +
  559 + // Return
  560 + return;
  561 +}
  562 +
  563 +
  564 +/***********************************************************************//**
  565 + * @brief Delete class members
  566 + ***************************************************************************/
  567 +void GModelTemporalFunc::free_members(void)
  568 +{
  569 + // Return
  570 + return;
  571 +}
  572 +
  573 +
  574 +/***********************************************************************//**
  575 + * @brief Load nodes from file
  576 + *
  577 + * @param[in] filename File name.
  578 + *
  579 + * @exception GException::invalid_value
  580 + * File function FITS file is invalid
  581 + *
  582 + * Load the temporal file function from a FITS file. The file function is
  583 + * expected in the first extension of the FITS file, containing two mandatory
  584 + * columns with names "TIME" and "NORM".
  585 + ***************************************************************************/
  586 +void GModelTemporalFunc::load_nodes(const GFilename& filename)
  587 +{
  588 + // Clear nodes and values
  589 + m_nodes.clear();
  590 + m_values.clear();
  591 +
  592 + // Set filename
  593 + m_filename = filename;
  594 +
  595 + // Load FITS file
  596 + GFits fits = GFits(filename);
  597 +
  598 + // Extract binary table (so far always load extension 1 as table)
  599 + GFitsTable* table = fits.table(1);
  600 +
  601 + // Read time reference from binary table
  602 + m_timeref.read(*table);
  603 +
  604 + // Extract columns
  605 + GFitsTableCol* time_col = (*table)["TIME"];
  606 + GFitsTableCol* time_norm = (*table)["NORM"];
  607 +
  608 + // Check that there are at least two nodes in table
  609 + if (time_col->length() < 2) {
  610 + std::string msg = "\"TIME\" column contains "+
  611 + gammalib::str(time_col->length())+" rows but at "
  612 + "least two rows are required. Please specify a valid "
  613 + "temporal file function.";
  614 + throw GException::invalid_value(G_LOAD_NODES, msg);
  615 + }
  616 +
  617 + // Check that both columns are consistent
  618 + if (time_col->length() != time_norm->length()) {
  619 + std::string msg = "\"TIME\" and \"NORM\" columns have inconsistent "
  620 + "number of rows ("+
  621 + gammalib::str(time_col->length())+", "+
  622 + gammalib::str(time_norm->length())+"). Please "
  623 + "specify a valid temporal file function.";
  624 + throw GException::invalid_value(G_LOAD_NODES, msg);
  625 + }
  626 +
  627 + // Extract nodes
  628 + for (int i = 0; i < time_col->length(); ++i) {
  629 + m_nodes.append(time_col->real(i));
  630 + m_values.push_back(time_norm->real(i));
  631 + }
  632 +
  633 + // Set minimum and maximum times (assumes that times are ordered)
  634 + m_tmin.set(time_col->real(0), m_timeref);
  635 + m_tmax.set(time_col->real(time_col->length()-1), m_timeref);
  636 +
  637 + // Close FITS file
  638 + fits.close();
  639 +
  640 + // Set pre-computation cache
  641 + //set_cache();
  642 +
  643 + // Return
  644 + return;
  645 +}
  646 +
  647 +
  648 +/***********************************************************************//**
  649 + * @brief Set MC pre-computation cache
  650 + *
  651 + * @param[in] tmin Minimum time of time interval.
  652 + * @param[in] tmax Maximum time of time interval.
  653 + *
  654 + * This method sets up an array of indices and the cumulative distribution
  655 + * function needed for MC simulations.
  656 + ***************************************************************************/
  657 +void GModelTemporalFunc::mc_update(const GTime& tmin,
  658 + const GTime& tmax) const
  659 +{
  660 + // Update the cache only if the time interval has changed
  661 + if (tmin != m_mc_tmin || tmax != m_mc_tmax) {
  662 +
  663 + // Store new time interval
  664 + m_mc_tmin = tmin;
  665 + m_mc_tmax = tmax;
  666 +
  667 + // Initialise cache
  668 + m_mc_rate = 0.0;
  669 + m_mc_cum.clear();
  670 + m_mc_slope.clear();
  671 + m_mc_offset.clear();
  672 + m_mc_time.clear();
  673 + m_mc_dt.clear();
  674 +
  675 + // Initialise duration
  676 + double duration = 0.0;
  677 +
  678 + // Continue only if time interval overlaps with temporal file function
  679 + // and if time interval is valid
  680 + if (tmax > m_tmin && tmin < m_tmax && tmax > tmin) {
  681 +
  682 + // Loop over all intervals between the nodes
  683 + for (int i = 0; i < m_nodes.size()-1; ++i) {
  684 +
  685 + // Get start and stop time of interval
  686 + GTime tstart(m_nodes[i], m_timeref);
  687 + GTime tstop(m_nodes[i+1], m_timeref);
  688 +
  689 + // Make sure that we only consider the interval that overlaps
  690 + // with the requested time interval
  691 + if (tmin > tstart) {
  692 + tstart = tmin;
  693 + }
  694 + if (tmax < tstop) {
  695 + tstop = tmax;
  696 + }
  697 +
  698 + // If time interval is empty then skip this node
  699 + if (tstart >= tstop) {
  700 + continue;
  701 + }
  702 +
  703 + // Evaluate rate at start and stop time
  704 + double rstart = eval(tstart);
  705 + double rstop = eval(tstop);
  706 +
  707 + // Compute mean normalization for this interval
  708 + double norm = 0.5 * (rstart + rstop);
  709 +
  710 + // Compute integral for this interval
  711 + double dt = tstop - tstart;
  712 + double cum = norm * dt;
  713 +
  714 + // Compute slope, offset and start time of interval
  715 + double slope = (rstop - rstart) / dt;
  716 + double offset = rstart;
  717 + double renorm = (0.5 * slope * dt + offset) * dt;
  718 + slope /= renorm;
  719 + offset /= renorm;
  720 +
  721 + // Update mean rate and duration
  722 + m_mc_rate += cum;
  723 + duration += dt;
  724 +
  725 + // Put mean normalization and integral on stack
  726 + m_mc_cum.push_back(cum);
  727 + m_mc_slope.push_back(slope);
  728 + m_mc_offset.push_back(offset);
  729 + m_mc_time.push_back(tstart.convert(m_timeref));
  730 + m_mc_dt.push_back(dt);
  731 +
  732 + } // endfor: next interval
  733 +
  734 + // Build cumulative distribution
  735 + for (int i = 1; i < m_mc_cum.size(); ++i) {
  736 + m_mc_cum[i] += m_mc_cum[i-1];
  737 + }
  738 + double norm = m_mc_cum[m_mc_cum.size()-1];
  739 + for (int i = 0; i < m_mc_cum.size(); ++i) {
  740 + m_mc_cum[i] /= norm;
  741 + }
  742 +
  743 + // Compute mean rate
  744 + if (duration > 0.0) {
  745 + m_mc_rate /= duration;
  746 + }
  747 +
  748 + } // endif: time interval overlaps
  749 +
  750 + } // endif: Update was required
  751 +
  752 + // Return
  753 + return;
  754 +}
... ...
src/model/Makefile.am
1 1 #############################################################################
2 2 # Makefile.am for model module #
3 3 # ------------------------------------------------------------------------- #
4   -# Copyright (C) 2016 Juergen Knoedlseder #
  4 +# Copyright (C) 2016-2017 Juergen Knoedlseder #
5 5 # ------------------------------------------------------------------------- #
6 6 # #
7 7 # This program is free software: you can redistribute it and/or modify #
... ... @@ -65,6 +65,7 @@ sources = GModelPar.cpp \
65 65 GModelTemporal.cpp \
66 66 GModelTemporalRegistry.cpp \
67 67 GModelTemporalConst.cpp \
  68 + GModelTemporalFunc.cpp \
68 69 GException_model.cpp
69 70  
70 71 # Build libtool library
... ...
test/data/model_temporal_filefct.xml 0 → 100644
  1 +<?xml version="1.0" standalone="no"?>
  2 +<source_library title="source library">
  3 + <source name="Crab" type="PointSource">
  4 + <spectrum type="PowerLaw">
  5 + <parameter name="Prefactor" scale="1e-16" value="5.7" min="1e-07" max="1000.0" free="1"/>
  6 + <parameter name="Index" scale="-1" value="2.48" min="0.0" max="+5.0" free="1"/>
  7 + <parameter name="PivotEnergy" scale="1e6" value="0.3" min="0.01" max="1000.0" free="0"/>
  8 + </spectrum>
  9 + <spatialModel type="PointSource">
  10 + <parameter free="0" max="360" min="-360" name="RA" scale="1" value="83.6331" />
  11 + <parameter free="0" max="90" min="-90" name="DEC" scale="1" value="22.0145" />
  12 + </spatialModel>
  13 + <temporal type="FileFunction" file="temp_filefunction.fits">
  14 + <parameter name="Normalization" scale="1" value="1.0" min="0.0" max="1000.0" free="0"/>
  15 + </temporal>
  16 + </source>
  17 +</source_library>
... ...
test/data/temp_filefunction.fits 0 → 100644
No preview for this file type
test/test_GModel.cpp
1 1 /***************************************************************************
2 2 * test_GModel.cpp - test GModel class *
3 3 * ----------------------------------------------------------------------- *
4   - * copyright (C) 2009-2016 by Juergen Knoedlseder *
  4 + * copyright (C) 2009-2017 by Juergen Knoedlseder *
5 5 * ----------------------------------------------------------------------- *
6 6 * *
7 7 * This program is free software: you can redistribute it and/or modify *
... ... @@ -48,10 +48,11 @@ void TestGModel::set(void)
48 48 std::string datadir = std::getenv("TEST_DATA");
49 49  
50 50 // Set test files
51   - m_map_file = datadir + "/cena_lobes_parkes.fits";
52   - m_cube_file = datadir + "/test_cube.fits";
53   - m_filefct = datadir + "/filefunction.txt";
54   - m_xml_file = datadir + "/crab.xml";
  51 + m_map_file = datadir + "/cena_lobes_parkes.fits";
  52 + m_cube_file = datadir + "/test_cube.fits";
  53 + m_filefct = datadir + "/filefunction.txt";
  54 + m_temp_filefct = datadir + "/temp_filefunction.fits";
  55 + m_xml_file = datadir + "/crab.xml";
55 56  
56 57 // Set model definiton XML files
57 58 m_xml_model_point_const = datadir + "/model_point_const.xml";
... ... @@ -77,6 +78,7 @@ void TestGModel::set(void)
77 78 m_xml_model_elliptical_disk = datadir + "/model_elliptical_disk.xml";
78 79 m_xml_model_elliptical_gauss = datadir + "/model_elliptical_gauss.xml";
79 80 m_xml_model_spatial_composite = datadir + "/model_spatial_composite.xml";
  81 + m_xml_model_point_temp_filefct = datadir + "/model_temporal_filefct.xml";
80 82  
81 83 // Set legacy model definition XML files
82 84 m_xml_legacy_radial_disk = datadir + "/legacy_radial_disk.xml";
... ... @@ -158,6 +160,8 @@ void TestGModel::set(void)
158 160 // Append temporal model tests
159 161 append(static_cast<pfunction>(&TestGModel::test_temp_const),
160 162 "Test GModelTemporalConst");
  163 + append(static_cast<pfunction>(&TestGModel::test_temp_filefct),
  164 + "Test GModelTemporalFunc");
161 165  
162 166 // Append genereric model tests
163 167 append(static_cast<pfunction>(&TestGModel::test_model), "Test GModel");
... ... @@ -2158,6 +2162,58 @@ void TestGModel::test_temp_const(void)
2158 2162  
2159 2163  
2160 2164 /***********************************************************************//**
  2165 + * @brief Test GModelTemporalFunc class
  2166 + ***************************************************************************/
  2167 +void TestGModel::test_temp_filefct(void)
  2168 +{
  2169 + // Test void constructor
  2170 + GModelTemporalFunc model1;
  2171 + test_value(model1.type(), "FileFunction", "Check void model type");
  2172 +
  2173 + // Test value constructor
  2174 + GModelTemporalFunc model2(m_temp_filefct, 2.0);
  2175 + test_value(model2.filename().url(), m_temp_filefct,
  2176 + "Check temporal file function data file name");
  2177 + test_value(model2.norm(), 2.0);
  2178 +
  2179 + // Test XML constructor
  2180 + GXml xml(m_xml_model_point_temp_filefct);
  2181 + GXmlElement* element = xml.element(0)->element(0)->element("temporal", 0);
  2182 + GModelTemporalFunc model3(*element);
  2183 + test_value(model3.size(), 1);
  2184 + test_value(model3.type(), "FileFunction", "Check model type");
  2185 + test_value(model3.filename().url(), m_temp_filefct,
  2186 + "Check temporal file function data file name");
  2187 + test_value(model3.norm(), 1.0);
  2188 +
  2189 + // Test filename method
  2190 + model3.filename(m_temp_filefct);
  2191 + test_value(model3.filename().url(), m_temp_filefct,
  2192 + "Check temporal file function data file name");
  2193 +
  2194 + // Test norm method
  2195 + model3.norm(3.0);
  2196 + test_value(model3.norm(), 3.0);
  2197 +
  2198 + // Test operator access
  2199 + const char* strarray[] = {"Normalization"};
  2200 + for (int i = 0; i < 1; ++i) {
  2201 + std::string keyname(strarray[i]);
  2202 + model3[keyname].remove_range(); // To allow setting of any value
  2203 + model3[keyname].value(2.1);
  2204 + model3[keyname].error(1.9);
  2205 + model3[keyname].gradient(0.8);
  2206 + test_value(model3[keyname].value(), 2.1);
  2207 + test_value(model3[keyname].error(), 1.9);
  2208 + test_value(model3[keyname].gradient(), 0.8);
  2209 + }
  2210 +
  2211 + // Return
  2212 + return;
  2213 +}
  2214 +
  2215 +
  2216 +/***********************************************************************//**
2161 2217 * @brief Test loading and saving of XML model.
2162 2218 *
2163 2219 * @param[in] name Model name.
... ...
test/test_GModel.hpp
1 1 /***************************************************************************
2 2 * test_GModel.hpp - Test model module *
3 3 * ----------------------------------------------------------------------- *
4   - * copyright (C) 2012-2016 by Jean-Baptiste Cayrou *
  4 + * copyright (C) 2012-2017 by Jean-Baptiste Cayrou *
5 5 * ----------------------------------------------------------------------- *
6 6 * *
7 7 * This program is free software: you can redistribute it and/or modify *
... ... @@ -76,6 +76,7 @@ public:
76 76 void test_spectral_composite(void);
77 77 void test_spectral_model(void);
78 78 void test_temp_const(void);
  79 + void test_temp_filefct(void);
79 80 void test_model(void);
80 81 void test_models(void);
81 82 void test_model_registry(void);
... ... @@ -101,6 +102,7 @@ private:
101 102 std::string m_map_file;
102 103 std::string m_cube_file;
103 104 std::string m_filefct;
  105 + std::string m_temp_filefct;
104 106 std::string m_xml_file;
105 107  
106 108 // Model definiton XML files
... ... @@ -127,6 +129,7 @@ private:
127 129 std::string m_xml_model_elliptical_disk;
128 130 std::string m_xml_model_elliptical_gauss;
129 131 std::string m_xml_model_spatial_composite;
  132 + std::string m_xml_model_point_temp_filefct;
130 133  
131 134 // Legacy model definition XML files
132 135 std::string m_xml_legacy_radial_disk;
... ...