From 5d30b7bb747d87a6931fc21ce99fe606c9d6f3a2 Mon Sep 17 00:00:00 2001 From: Jürgen Knödlseder Date: Wed, 25 Jan 2017 17:18:29 +0100 Subject: [PATCH] Add GModelTemporalFunc class (#1914) --- ChangeLog | 3 ++- NEWS | 8 +++++++- include/GModelTemporalFunc.hpp | 207 +++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++ include/GammaLib.hpp | 3 ++- include/Makefile.am | 5 +++-- pyext/GModelTemporalFunc.i | 75 +++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++ pyext/gammalib/model.i | 3 ++- src/model/GModelTemporalFunc.cpp | 754 ++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++ src/model/Makefile.am | 3 ++- test/data/model_temporal_filefct.xml | 17 +++++++++++++++++ test/data/temp_filefunction.fits | Bin 0 -> 8640 bytes test/test_GModel.cpp | 66 +++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++----- test/test_GModel.hpp | 5 ++++- 13 files changed, 1136 insertions(+), 13 deletions(-) create mode 100644 include/GModelTemporalFunc.hpp create mode 100644 pyext/GModelTemporalFunc.i create mode 100644 src/model/GModelTemporalFunc.cpp create mode 100644 test/data/model_temporal_filefct.xml create mode 100644 test/data/temp_filefunction.fits diff --git a/ChangeLog b/ChangeLog index 330d67f..21f11e8 100644 --- a/ChangeLog +++ b/ChangeLog @@ -1,8 +1,9 @@ -2017-01-24 +2017-01-25 * Version 1.2.0 released ========================= + Add GModelTemporalFunc call (#1914) Add GSkyMap FITS HDU constructor Add leap second for 1-1-2017 (#1887) Add region() method to GModelSpatial (#1764) diff --git a/NEWS b/NEWS index 377f4e5..8fb09ab 100644 --- a/NEWS +++ b/NEWS @@ -1,6 +1,6 @@ New Features and Important Changes in GammaLib 1.2.0 -24 January 2017 +25 January 2017 1. Introduction @@ -17,6 +17,7 @@ The following classes have been added: - GModelSpectralMultiplicative - GModelSpectralComposite - GModelSpatialComposite +- GModelTemporalFunc - GVOTable - GNdarray - GFft @@ -169,6 +170,11 @@ derived classes that returns a pointer to a sky region that bounds the spatial model. For the moment, only circular sky regions have been implemented (#1764). +The GModelTemporalFunc class has been added to implement a temporal file +function for the simulation of transient sources. A unit test has been +added and the mc() and eval() methods have been validated using a dedicated +test script (#1914). + 9. Numerics module ------------------ diff --git a/include/GModelTemporalFunc.hpp b/include/GModelTemporalFunc.hpp new file mode 100644 index 0000000..cfa2bef --- /dev/null +++ b/include/GModelTemporalFunc.hpp @@ -0,0 +1,207 @@ +/*************************************************************************** + * GModelTemporalFunc.hpp - Temporal file function model class * + * ----------------------------------------------------------------------- * + * copyright (C) 2017 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 GModelTemporalFunc.hpp + * @brief File function temporal model class interface definition + * @author Juergen Knoedlseder + */ + +#ifndef GMODELTEMPORALFUNC_HPP +#define GMODELTEMPORALFUNC_HPP + +/* __ Includes ___________________________________________________________ */ +#include +#include +#include "GModelPar.hpp" +#include "GModelTemporal.hpp" +#include "GNodeArray.hpp" +#include "GFilename.hpp" +#include "GTimeReference.hpp" +#include "GTime.hpp" + +/* __ Forward declarations _______________________________________________ */ +class GRan; +class GXmlElement; + + +/***********************************************************************//** + * @class GModelTemporalFunc + * + * @brief File function temporal model class + * + * This class implements a light curve defined by nodes given specified in + * a FITS file. The model is defined by + * + * \f[ + * S_{\rm t}(t) = r(t) \times {\tt m\_norm} + * \f] + * + * where + * \f$r(t)$\f is the rate defined by linear interpolation between the nodes + * in a FITS file, and + * \f${\tt m\_norm}$\f is a normalisation constant. + ***************************************************************************/ +class GModelTemporalFunc : public GModelTemporal { + +public: + // Constructors and destructors + GModelTemporalFunc(void); + explicit GModelTemporalFunc(const GXmlElement& xml); + GModelTemporalFunc(const GFilename& filename, + const double& norm = 1.0); + GModelTemporalFunc(const GModelTemporalFunc& model); + virtual ~GModelTemporalFunc(void); + + // Operators + virtual GModelTemporalFunc& operator=(const GModelTemporalFunc& model); + + // Implemented virtual base class methods + virtual void clear(void); + virtual GModelTemporalFunc* clone(void) const; + virtual std::string classname(void) const; + virtual std::string type(void) const; + virtual double eval(const GTime& srcTime, + const bool& gradients = false) const; + virtual GTimes mc(const double& rate, const GTime& tmin, + const GTime& tmax, GRan& ran) const; + virtual void read(const GXmlElement& xml); + virtual void write(GXmlElement& xml) const; + virtual std::string print(const GChatter& chatter = NORMAL) const; + + // Other methods + const GFilename& filename(void) const; + void filename(const GFilename& filename); + double norm(void) const; + void norm(const double& norm); + +protected: + // Protected methods + void init_members(void); + void copy_members(const GModelTemporalFunc& model); + void free_members(void); + void load_nodes(const GFilename& filename); + void mc_update(const GTime& tmin, const GTime& tmax) const; + + // Protected members + GModelPar m_norm; //!< Normalization factor + mutable GNodeArray m_nodes; //!< Time nodes of function + std::vector m_values; //!< Function values at nodes + GFilename m_filename; //!< Name of file function + GTimeReference m_timeref; //!< Time reference + GTime m_tmin; //!< Minimum time of model + GTime m_tmax; //!< Maximum time of model + + // Cached members for MC + mutable GTime m_mc_tmin; //!< Minimum time + mutable GTime m_mc_tmax; //!< Maximum time + mutable double m_mc_rate; //!< Mean rate + mutable std::vector m_mc_cum; //!< Cumulative distribution + mutable std::vector m_mc_slope; //!< Slope of interval + mutable std::vector m_mc_offset; //!< Offset of interval + mutable std::vector m_mc_time; //!< Start time of interval + mutable std::vector m_mc_dt; //!< Length of interval +}; + + +/***********************************************************************//** + * @brief Return class name + * + * @return String containing the class name ("GModelTemporalFunc"). + ***************************************************************************/ +inline +std::string GModelTemporalFunc::classname(void) const +{ + return ("GModelTemporalFunc"); +} + + +/***********************************************************************//** + * @brief Return model type + * + * @return "FileFunction". + * + * Returns the type of the temporal model. + ***************************************************************************/ +inline +std::string GModelTemporalFunc::type(void) const +{ + return "FileFunction"; +} + + +/***********************************************************************//** + * @brief Return normalization factor + * + * @return Normalization factor. + * + * Returns the normalization factor. + ***************************************************************************/ +inline +double GModelTemporalFunc::norm(void) const +{ + return (m_norm.value()); +} + + +/***********************************************************************//** + * @brief Set normalization factor + * + * @param[in] norm Normalization factor. + * + * Sets the normalization factor. + ***************************************************************************/ +inline +void GModelTemporalFunc::norm(const double& norm) +{ + m_norm.value(norm); + return; +} + + +/***********************************************************************//** + * @brief Return file name + * + * @return Name of node file. + * + * Returns the name of the file function node file. + ***************************************************************************/ +inline +const GFilename& GModelTemporalFunc::filename(void) const +{ + return (m_filename); +} + + +/***********************************************************************//** + * @brief Loads nodes from node file and set filename + * + * @param[in] filename Node file name. + * + * Loads the nodes from a file function node file and sets the filename. + ***************************************************************************/ +inline +void GModelTemporalFunc::filename(const GFilename& filename) +{ + load_nodes(filename); + return; +} + +#endif /* GMODELTEMPORALFUNC_HPP */ diff --git a/include/GammaLib.hpp b/include/GammaLib.hpp index 4e066e4..ceda7d7 100644 --- a/include/GammaLib.hpp +++ b/include/GammaLib.hpp @@ -1,7 +1,7 @@ /*************************************************************************** * GammaLib.hpp - Gamma-Ray Astronomy Library Header file * * ----------------------------------------------------------------------- * - * copyright (C) 2006-2016 by Juergen Knoedlseder * + * copyright (C) 2006-2017 by Juergen Knoedlseder * * ----------------------------------------------------------------------- * * * * This program is free software: you can redistribute it and/or modify * @@ -228,5 +228,6 @@ #include "GModelTemporal.hpp" #include "GModelTemporalRegistry.hpp" #include "GModelTemporalConst.hpp" +#include "GModelTemporalFunc.hpp" #endif /* GAMMALIB_HPP */ diff --git a/include/Makefile.am b/include/Makefile.am index e8a291f..770a445 100644 --- a/include/Makefile.am +++ b/include/Makefile.am @@ -1,7 +1,7 @@ ############################################################################# # Makefile.am for include module # # ------------------------------------------------------------------------- # -# Copyright (C) 2016 Juergen Knoedlseder # +# Copyright (C) 2016-2017 Juergen Knoedlseder # # ------------------------------------------------------------------------- # # # # This program is free software: you can redistribute it and/or modify # @@ -185,4 +185,5 @@ pkginclude_HEADERS = GammaLib.hpp \ GModelSpectralComposite.hpp \ GModelTemporal.hpp \ GModelTemporalRegistry.hpp \ - GModelTemporalConst.hpp + GModelTemporalConst.hpp \ + GModelTemporalFunc.hpp diff --git a/pyext/GModelTemporalFunc.i b/pyext/GModelTemporalFunc.i new file mode 100644 index 0000000..ce07aec --- /dev/null +++ b/pyext/GModelTemporalFunc.i @@ -0,0 +1,75 @@ +/*************************************************************************** + * GModelTemporalFunc.i - Temporal file function model class * + * ----------------------------------------------------------------------- * + * copyright (C) 2017 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 GModelTemporalFunc.i + * @brief File function temporal model class interface definition + * @author Juergen Knoedlseder + */ +%{ +/* Put headers and other declarations here that are needed for compilation */ +#include "GModelTemporalFunc.hpp" +%} + + +/***********************************************************************//** + * @class GModelTemporalFunc + * + * @brief File function temporal model class + ***************************************************************************/ +class GModelTemporalFunc : public GModelTemporal { + +public: + // Constructors and destructors + GModelTemporalFunc(void); + explicit GModelTemporalFunc(const GXmlElement& xml); + GModelTemporalFunc(const GFilename& filename, + const double& norm = 1.0); + GModelTemporalFunc(const GModelTemporalFunc& model); + virtual ~GModelTemporalFunc(void); + + // Implemented virtual base class methods + virtual void clear(void); + virtual GModelTemporalFunc* clone(void) const; + virtual std::string classname(void) const; + virtual std::string type(void) const; + virtual double eval(const GTime& srcTime, + const bool& gradients = false) const; + virtual GTimes mc(const double& rate, const GTime& tmin, + const GTime& tmax, GRan& ran) const; + virtual void read(const GXmlElement& xml); + virtual void write(GXmlElement& xml) const; + + // Other methods + const GFilename& filename(void) const; + void filename(const GFilename& filename); + double norm(void) const; + void norm(const double& norm); +}; + + +/***********************************************************************//** + * @brief GModelTemporalFunc class extension + ***************************************************************************/ +%extend GModelTemporalFunc { + GModelTemporalFunc copy() { + return (*self); + } +}; diff --git a/pyext/gammalib/model.i b/pyext/gammalib/model.i index 003e6c3..67bba3c 100644 --- a/pyext/gammalib/model.i +++ b/pyext/gammalib/model.i @@ -1,7 +1,7 @@ /*************************************************************************** * model.i - Model module * * ----------------------------------------------------------------------- * - * copyright (C) 2011-2016 by Juergen Knoedlseder * + * copyright (C) 2011-2017 by Juergen Knoedlseder * * ----------------------------------------------------------------------- * * * * This program is free software: you can redistribute it and/or modify * @@ -160,3 +160,4 @@ %include "GModelTemporal.i" %include "GModelTemporalRegistry.i" %include "GModelTemporalConst.i" +%include "GModelTemporalFunc.i" diff --git a/src/model/GModelTemporalFunc.cpp b/src/model/GModelTemporalFunc.cpp new file mode 100644 index 0000000..f7340d0 --- /dev/null +++ b/src/model/GModelTemporalFunc.cpp @@ -0,0 +1,754 @@ +/*************************************************************************** + * GModelTemporalFunc.cpp - Temporal file function model class * + * ----------------------------------------------------------------------- * + * copyright (C) 2017 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 GModelTemporalFunc.cpp + * @brief File function temporal model class interface implementation + * @author Juergen Knoedlseder + */ + +/* __ Includes ___________________________________________________________ */ +#ifdef HAVE_CONFIG_H +#include +#endif +#include "GTools.hpp" +#include "GException.hpp" +#include "GModelTemporalFunc.hpp" +#include "GModelTemporalRegistry.hpp" +#include "GFits.hpp" + +/* __ Constants __________________________________________________________ */ + +/* __ Globals ____________________________________________________________ */ +const GModelTemporalFunc g_temporal_func_seed; +const GModelTemporalRegistry g_temporal_func_registry(&g_temporal_func_seed); + +/* __ Method name definitions ____________________________________________ */ +#define G_READ "GModelTemporalFunc::read(GXmlElement&)" +#define G_WRITE "GModelTemporalFunc::write(GXmlElement&)" +#define G_LOAD_NODES "GModelTemporalFunc::load_nodes(GFilename&)" + +/* __ Macros _____________________________________________________________ */ + +/* __ Coding definitions _________________________________________________ */ + +/* __ Debug definitions __________________________________________________ */ + + +/*========================================================================== + = = + = Constructors/destructors = + = = + ==========================================================================*/ + +/***********************************************************************//** + * @brief Void constructor + ***************************************************************************/ +GModelTemporalFunc::GModelTemporalFunc(void) : GModelTemporal() +{ + // Initialise members + init_members(); + + // Return + return; +} + + +/***********************************************************************//** + * @brief File constructor + * + * @param[in] filename File name of nodes. + * @param[in] norm Normalization factor. + * + * Constructs spectral file function model from a list of nodes that is found + * in the specified FITS file. See the load_nodes() method for more + * information about the expected structure of the file. + ***************************************************************************/ +GModelTemporalFunc::GModelTemporalFunc(const GFilename& filename, + const double& norm) : + GModelTemporal() +{ + // Initialise members + init_members(); + + // Load nodes + load_nodes(filename); + + // Set normalization + m_norm.value(norm); + + // Return + return; +} + + +/***********************************************************************//** + * @brief XML constructor + * + * @param[in] xml XML element. + * + * Constructs constant temporal model by extracting information from an XML + * element. See the read() method for more information about the expected + * structure of the XML element. + ***************************************************************************/ +GModelTemporalFunc::GModelTemporalFunc(const GXmlElement& xml) : + GModelTemporal() +{ + // Initialise members + init_members(); + + // Read information from XML element + read(xml); + + // Return + return; +} + + +/***********************************************************************//** + * @brief Copy constructor + * + * @param[in] model File function temporal model. + ***************************************************************************/ +GModelTemporalFunc::GModelTemporalFunc(const GModelTemporalFunc& model) : + GModelTemporal(model) +{ + // Initialise members + init_members(); + + // Copy members + copy_members(model); + + // Return + return; +} + + +/***********************************************************************//** + * @brief Destructor + ***************************************************************************/ +GModelTemporalFunc::~GModelTemporalFunc(void) +{ + // Free members + free_members(); + + // Return + return; +} + + +/*========================================================================== + = = + = Operators = + = = + ==========================================================================*/ + +/***********************************************************************//** + * @brief Assignment operator + * + * @param[in] model File function temporal model. + * @return File function temporal model. + ***************************************************************************/ +GModelTemporalFunc& GModelTemporalFunc::operator=(const GModelTemporalFunc& model) +{ + // Execute only if object is not identical + if (this != &model) { + + // Copy base class members + this->GModelTemporal::operator=(model); + + // Free members + free_members(); + + // Initialise members + init_members(); + + // Copy members + copy_members(model); + + } // endif: object was not identical + + // Return + return *this; +} + + +/*========================================================================== + = = + = Public methods = + = = + ==========================================================================*/ + +/***********************************************************************//** + * @brief Clear file function temporal model + ***************************************************************************/ +void GModelTemporalFunc::clear(void) +{ + // Free class members (base and derived classes, derived class first) + free_members(); + this->GModelTemporal::free_members(); + + // Initialise members + this->GModelTemporal::init_members(); + init_members(); + + // Return + return; +} + + +/***********************************************************************//** + * @brief Clone file function temporal model + * + * @return Pointer to deep copy of file function temporal model. + ***************************************************************************/ +GModelTemporalFunc* GModelTemporalFunc::clone(void) const +{ + // Clone constant temporal model + return new GModelTemporalFunc(*this); +} + + +/***********************************************************************//** + * @brief Evaluate function + * + * @param[in] srcTime True photon arrival time (not used). + * @param[in] gradients Compute gradients? + * @return Value of temporal model. + * + * Computes + * + * \f[ + * S_{\rm t}(t) = r(t) \times {\tt m\_norm} + * \f] + * + * where + * \f$r(t)\f$ is the linear interpolation between the nodes of the file + * function and + * \f${\tt m\_norm}\f$ is the normalization constant. + * + * If the @p gradients flag is true the method will also evaluate the partial + * derivatives of the model with respect to the normalization parameter using + * + * \f[ + * \frac{\delta S_{\rm t}(t)}{\delta {\tt m\_norm}} = r(t) + * \f] + ***************************************************************************/ +double GModelTemporalFunc::eval(const GTime& srcTime, + const bool& gradients) const +{ + // Initialise value + double value = 0.0; + + // Optionally initialise gradient + if (gradients) { + m_norm.factor_gradient(0.0); + } + + // Set value if the time is within the validity range + if (srcTime >= m_tmin && srcTime <= m_tmax) { + + // Interpolate file function + double func = m_nodes.interpolate(srcTime.convert(m_timeref), m_values); + + // Compute function value + value = m_norm.value() * func; + + // Optionally compute gradients + if (gradients) { + + // Compute partial derivatives of the parameter values + double g_norm = (m_norm.is_free()) ? m_norm.scale() * func : 0.0; + + // Set gradients + m_norm.factor_gradient(g_norm); + + } // endif: gradient computation was requested + + } // endif: time was within the validity range + + // Return + return value; +} + + +/***********************************************************************//** + * @brief Returns vector of random event times + * + * @param[in] rate Mean event rate (events per second). + * @param[in] tmin Minimum event time. + * @param[in] tmax Maximum event time. + * @param[in,out] ran Random number generator. + * + * This method returns a vector of random event times between @p tmin and + * @p tmax assuming a light curve specified by a file function. + ***************************************************************************/ +GTimes GModelTemporalFunc::mc(const double& rate, const GTime& tmin, + const GTime& tmax, GRan& ran) const +{ + // Allocates empty vector of times + GTimes times; + + // Update Monte Carlo cache + mc_update(tmin, tmax); + + // Continue only if mean rate is positive and cache is not empty + if (m_mc_rate > 0.0 && m_mc_cum.size() > 0) { + + // Compute mean event rate (in events per seconds) + double lambda = rate * norm() * m_mc_rate; + + // Compute number of times to be sampled + int ntimes = int(ran.poisson(lambda)+0.5); + + // Loop over number of times + for (int i = 0; i < ntimes; ++i) { + + // Determine in which bin we reside + int inx = 0; + if (m_mc_cum.size() > 1) { + double u = ran.uniform(); + for (inx = m_mc_cum.size()-1; inx > 0; --inx) { + if (m_mc_cum[inx-1] <= u) { + break; + } + } + } + + // Get random time + double seconds; + if (m_mc_slope[inx] == 0.0) { + seconds = m_mc_dt[inx] * ran.uniform() + m_mc_time[inx]; + } + else { + seconds = (std::sqrt(m_mc_offset[inx]*m_mc_offset[inx] + + 2.0 * m_mc_slope[inx] * ran.uniform()) - + m_mc_offset[inx]) / m_mc_slope[inx] + m_mc_time[inx]; + } + + // Append random time + GTime time(seconds, m_timeref); + times.append(time); + + } // endfor: looped over times + + } // endif: cache was valid + + // Return vector of times + return times; +} + + +/***********************************************************************//** + * @brief Read model from XML element + * + * @param[in] xml XML element. + * + * Writes the temporal information from an XML element having the format + * + * + * + * + ***************************************************************************/ +void GModelTemporalFunc::read(const GXmlElement& xml) +{ + // Get parameter pointers + const GXmlElement* norm = gammalib::xml_get_par(G_READ, xml, m_norm.name()); + + // Read parameters + m_norm.read(*norm); + + // Load nodes from file + load_nodes(gammalib::xml_file_expand(xml, xml.attribute("file"))); + + // Return + return; +} + + +/***********************************************************************//** + * @brief Write model into XML element + * + * @param[in] xml XML element. + * + * @exception GExpection::invalid_value + * Invalid XML format encountered. + * + * Writes the temporal information into an XML element in the format + * + * + * + * + ***************************************************************************/ +void GModelTemporalFunc::write(GXmlElement& xml) const +{ + // Set model type + if (xml.attribute("type") == "") { + xml.attribute("type", type()); + } + + // Verify model type + if (xml.attribute("type") != type()) { + std::string msg = "Temporal model of type "+xml.attribute("type")+ + " encountered while \""+type()+"\" was expected."; + throw GException::invalid_value(G_WRITE, msg); + } + + // Get XML parameters + GXmlElement* norm = gammalib::xml_need_par(G_WRITE, xml, m_norm.name()); + + // Write parameters + m_norm.write(*norm); + + // Set file attribute + xml.attribute("file", gammalib::xml_file_reduce(xml, m_filename)); + + // Return + return; +} + + +/***********************************************************************//** + * @brief Print file function information + * + * @param[in] chatter Chattiness. + * @return String containing model information. + ***************************************************************************/ +std::string GModelTemporalFunc::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("=== GModelTemporalFunc ==="); + + // Append information + result.append("\n"+gammalib::parformat("Function file")); + result.append(m_filename.url()); + result.append("\n"+gammalib::parformat("Number of parameters")); + result.append(gammalib::str(size())); + for (int i = 0; i < size(); ++i) { + result.append("\n"+m_pars[i]->print(chatter)); + } + + // Append reference MJD and time range + result.append("\n"+gammalib::parformat("Reference MDJ")); + result.append(gammalib::str(m_timeref.mjdref())); + result.append("\n"+gammalib::parformat("Time range")); + result.append(gammalib::str(m_tmin.convert(m_timeref))); + result.append(" - "); + result.append(gammalib::str(m_tmax.convert(m_timeref))); + result.append(" "+m_timeref.timeunit()); + result.append(" ("+m_timeref.timesys()+")"); + + // Append node information + result.append("\n"+gammalib::parformat("Number of nodes")); + result.append(gammalib::str(m_nodes.size())); + for (int i = 0; i < m_nodes.size(); ++i) { + result.append("\n"); + result.append(gammalib::parformat(" Time "+gammalib::str(m_nodes[i])+" s")); + result.append(gammalib::str(m_values[i])); + } + + } // endif: chatter was not silent + + // Return result + return result; +} + + +/*========================================================================== + = = + = Private methods = + = = + ==========================================================================*/ + +/***********************************************************************//** + * @brief Initialise class members + ***************************************************************************/ +void GModelTemporalFunc::init_members(void) +{ + // Initialise normalisation parameter + m_norm.clear(); + m_norm.name("Normalization"); + m_norm.unit("(relative value)"); + m_norm.scale(1.0); + m_norm.value(1.0); + m_norm.range(0.0,1000.0); + m_norm.fix(); + m_norm.gradient(0.0); + m_norm.has_grad(true); + + // Set parameter pointer(s) + m_pars.clear(); + m_pars.push_back(&m_norm); + + // Initialise members + m_nodes.clear(); + m_values.clear(); + m_filename.clear(); + m_timeref.clear(); + m_tmin.clear(); + m_tmax.clear(); + + // Initialise cache + m_mc_tmin.clear(); + m_mc_tmax.clear(); + m_mc_rate = 0.0; + m_mc_cum.clear(); + m_mc_slope.clear(); + m_mc_offset.clear(); + m_mc_time.clear(); + m_mc_dt.clear(); + + // Return + return; +} + + +/***********************************************************************//** + * @brief Copy class members + * + * @param[in] model Constant temporal model + ***************************************************************************/ +void GModelTemporalFunc::copy_members(const GModelTemporalFunc& model) +{ + // Copy members + m_norm = model.m_norm; + m_nodes = model.m_nodes; + m_values = model.m_values; + m_filename = model.m_filename; + m_timeref = model.m_timeref; + m_tmin = model.m_tmin; + m_tmax = model.m_tmax; + + // Copy cache + m_mc_tmin = model.m_mc_tmin; + m_mc_tmax = model.m_mc_tmax; + m_mc_rate = model.m_mc_rate; + m_mc_cum = model.m_mc_cum; + m_mc_slope = model.m_mc_slope; + m_mc_offset = model.m_mc_offset; + m_mc_time = model.m_mc_time; + m_mc_dt = model.m_mc_dt; + + // Set parameter pointer(s) + m_pars.clear(); + m_pars.push_back(&m_norm); + + // Return + return; +} + + +/***********************************************************************//** + * @brief Delete class members + ***************************************************************************/ +void GModelTemporalFunc::free_members(void) +{ + // Return + return; +} + + +/***********************************************************************//** + * @brief Load nodes from file + * + * @param[in] filename File name. + * + * @exception GException::invalid_value + * File function FITS file is invalid + * + * Load the temporal file function from a FITS file. The file function is + * expected in the first extension of the FITS file, containing two mandatory + * columns with names "TIME" and "NORM". + ***************************************************************************/ +void GModelTemporalFunc::load_nodes(const GFilename& filename) +{ + // Clear nodes and values + m_nodes.clear(); + m_values.clear(); + + // Set filename + m_filename = filename; + + // Load FITS file + GFits fits = GFits(filename); + + // Extract binary table (so far always load extension 1 as table) + GFitsTable* table = fits.table(1); + + // Read time reference from binary table + m_timeref.read(*table); + + // Extract columns + GFitsTableCol* time_col = (*table)["TIME"]; + GFitsTableCol* time_norm = (*table)["NORM"]; + + // Check that there are at least two nodes in table + if (time_col->length() < 2) { + std::string msg = "\"TIME\" column contains "+ + gammalib::str(time_col->length())+" rows but at " + "least two rows are required. Please specify a valid " + "temporal file function."; + throw GException::invalid_value(G_LOAD_NODES, msg); + } + + // Check that both columns are consistent + if (time_col->length() != time_norm->length()) { + std::string msg = "\"TIME\" and \"NORM\" columns have inconsistent " + "number of rows ("+ + gammalib::str(time_col->length())+", "+ + gammalib::str(time_norm->length())+"). Please " + "specify a valid temporal file function."; + throw GException::invalid_value(G_LOAD_NODES, msg); + } + + // Extract nodes + for (int i = 0; i < time_col->length(); ++i) { + m_nodes.append(time_col->real(i)); + m_values.push_back(time_norm->real(i)); + } + + // Set minimum and maximum times (assumes that times are ordered) + m_tmin.set(time_col->real(0), m_timeref); + m_tmax.set(time_col->real(time_col->length()-1), m_timeref); + + // Close FITS file + fits.close(); + + // Set pre-computation cache + //set_cache(); + + // Return + return; +} + + +/***********************************************************************//** + * @brief Set MC pre-computation cache + * + * @param[in] tmin Minimum time of time interval. + * @param[in] tmax Maximum time of time interval. + * + * This method sets up an array of indices and the cumulative distribution + * function needed for MC simulations. + ***************************************************************************/ +void GModelTemporalFunc::mc_update(const GTime& tmin, + const GTime& tmax) const +{ + // Update the cache only if the time interval has changed + if (tmin != m_mc_tmin || tmax != m_mc_tmax) { + + // Store new time interval + m_mc_tmin = tmin; + m_mc_tmax = tmax; + + // Initialise cache + m_mc_rate = 0.0; + m_mc_cum.clear(); + m_mc_slope.clear(); + m_mc_offset.clear(); + m_mc_time.clear(); + m_mc_dt.clear(); + + // Initialise duration + double duration = 0.0; + + // Continue only if time interval overlaps with temporal file function + // and if time interval is valid + if (tmax > m_tmin && tmin < m_tmax && tmax > tmin) { + + // Loop over all intervals between the nodes + for (int i = 0; i < m_nodes.size()-1; ++i) { + + // Get start and stop time of interval + GTime tstart(m_nodes[i], m_timeref); + GTime tstop(m_nodes[i+1], m_timeref); + + // Make sure that we only consider the interval that overlaps + // with the requested time interval + if (tmin > tstart) { + tstart = tmin; + } + if (tmax < tstop) { + tstop = tmax; + } + + // If time interval is empty then skip this node + if (tstart >= tstop) { + continue; + } + + // Evaluate rate at start and stop time + double rstart = eval(tstart); + double rstop = eval(tstop); + + // Compute mean normalization for this interval + double norm = 0.5 * (rstart + rstop); + + // Compute integral for this interval + double dt = tstop - tstart; + double cum = norm * dt; + + // Compute slope, offset and start time of interval + double slope = (rstop - rstart) / dt; + double offset = rstart; + double renorm = (0.5 * slope * dt + offset) * dt; + slope /= renorm; + offset /= renorm; + + // Update mean rate and duration + m_mc_rate += cum; + duration += dt; + + // Put mean normalization and integral on stack + m_mc_cum.push_back(cum); + m_mc_slope.push_back(slope); + m_mc_offset.push_back(offset); + m_mc_time.push_back(tstart.convert(m_timeref)); + m_mc_dt.push_back(dt); + + } // endfor: next interval + + // Build cumulative distribution + for (int i = 1; i < m_mc_cum.size(); ++i) { + m_mc_cum[i] += m_mc_cum[i-1]; + } + double norm = m_mc_cum[m_mc_cum.size()-1]; + for (int i = 0; i < m_mc_cum.size(); ++i) { + m_mc_cum[i] /= norm; + } + + // Compute mean rate + if (duration > 0.0) { + m_mc_rate /= duration; + } + + } // endif: time interval overlaps + + } // endif: Update was required + + // Return + return; +} diff --git a/src/model/Makefile.am b/src/model/Makefile.am index 852827e..199f5e0 100644 --- a/src/model/Makefile.am +++ b/src/model/Makefile.am @@ -1,7 +1,7 @@ ############################################################################# # Makefile.am for model module # # ------------------------------------------------------------------------- # -# Copyright (C) 2016 Juergen Knoedlseder # +# Copyright (C) 2016-2017 Juergen Knoedlseder # # ------------------------------------------------------------------------- # # # # This program is free software: you can redistribute it and/or modify # @@ -65,6 +65,7 @@ sources = GModelPar.cpp \ GModelTemporal.cpp \ GModelTemporalRegistry.cpp \ GModelTemporalConst.cpp \ + GModelTemporalFunc.cpp \ GException_model.cpp # Build libtool library diff --git a/test/data/model_temporal_filefct.xml b/test/data/model_temporal_filefct.xml new file mode 100644 index 0000000..f9a4aa7 --- /dev/null +++ b/test/data/model_temporal_filefct.xml @@ -0,0 +1,17 @@ + + + + + + + + + + + + + + + + + diff --git a/test/data/temp_filefunction.fits b/test/data/temp_filefunction.fits new file mode 100644 index 0000000..9328017 Binary files /dev/null and b/test/data/temp_filefunction.fits differ diff --git a/test/test_GModel.cpp b/test/test_GModel.cpp index e60eb71..b9e99f3 100644 --- a/test/test_GModel.cpp +++ b/test/test_GModel.cpp @@ -1,7 +1,7 @@ /*************************************************************************** * test_GModel.cpp - test GModel class * * ----------------------------------------------------------------------- * - * copyright (C) 2009-2016 by Juergen Knoedlseder * + * copyright (C) 2009-2017 by Juergen Knoedlseder * * ----------------------------------------------------------------------- * * * * This program is free software: you can redistribute it and/or modify * @@ -48,10 +48,11 @@ void TestGModel::set(void) std::string datadir = std::getenv("TEST_DATA"); // Set test files - m_map_file = datadir + "/cena_lobes_parkes.fits"; - m_cube_file = datadir + "/test_cube.fits"; - m_filefct = datadir + "/filefunction.txt"; - m_xml_file = datadir + "/crab.xml"; + m_map_file = datadir + "/cena_lobes_parkes.fits"; + m_cube_file = datadir + "/test_cube.fits"; + m_filefct = datadir + "/filefunction.txt"; + m_temp_filefct = datadir + "/temp_filefunction.fits"; + m_xml_file = datadir + "/crab.xml"; // Set model definiton XML files m_xml_model_point_const = datadir + "/model_point_const.xml"; @@ -77,6 +78,7 @@ void TestGModel::set(void) m_xml_model_elliptical_disk = datadir + "/model_elliptical_disk.xml"; m_xml_model_elliptical_gauss = datadir + "/model_elliptical_gauss.xml"; m_xml_model_spatial_composite = datadir + "/model_spatial_composite.xml"; + m_xml_model_point_temp_filefct = datadir + "/model_temporal_filefct.xml"; // Set legacy model definition XML files m_xml_legacy_radial_disk = datadir + "/legacy_radial_disk.xml"; @@ -158,6 +160,8 @@ void TestGModel::set(void) // Append temporal model tests append(static_cast(&TestGModel::test_temp_const), "Test GModelTemporalConst"); + append(static_cast(&TestGModel::test_temp_filefct), + "Test GModelTemporalFunc"); // Append genereric model tests append(static_cast(&TestGModel::test_model), "Test GModel"); @@ -2158,6 +2162,58 @@ void TestGModel::test_temp_const(void) /***********************************************************************//** + * @brief Test GModelTemporalFunc class + ***************************************************************************/ +void TestGModel::test_temp_filefct(void) +{ + // Test void constructor + GModelTemporalFunc model1; + test_value(model1.type(), "FileFunction", "Check void model type"); + + // Test value constructor + GModelTemporalFunc model2(m_temp_filefct, 2.0); + test_value(model2.filename().url(), m_temp_filefct, + "Check temporal file function data file name"); + test_value(model2.norm(), 2.0); + + // Test XML constructor + GXml xml(m_xml_model_point_temp_filefct); + GXmlElement* element = xml.element(0)->element(0)->element("temporal", 0); + GModelTemporalFunc model3(*element); + test_value(model3.size(), 1); + test_value(model3.type(), "FileFunction", "Check model type"); + test_value(model3.filename().url(), m_temp_filefct, + "Check temporal file function data file name"); + test_value(model3.norm(), 1.0); + + // Test filename method + model3.filename(m_temp_filefct); + test_value(model3.filename().url(), m_temp_filefct, + "Check temporal file function data file name"); + + // Test norm method + model3.norm(3.0); + test_value(model3.norm(), 3.0); + + // Test operator access + const char* strarray[] = {"Normalization"}; + for (int i = 0; i < 1; ++i) { + std::string keyname(strarray[i]); + model3[keyname].remove_range(); // To allow setting of any value + model3[keyname].value(2.1); + model3[keyname].error(1.9); + model3[keyname].gradient(0.8); + test_value(model3[keyname].value(), 2.1); + test_value(model3[keyname].error(), 1.9); + test_value(model3[keyname].gradient(), 0.8); + } + + // Return + return; +} + + +/***********************************************************************//** * @brief Test loading and saving of XML model. * * @param[in] name Model name. diff --git a/test/test_GModel.hpp b/test/test_GModel.hpp index 9a2a693..1683626 100644 --- a/test/test_GModel.hpp +++ b/test/test_GModel.hpp @@ -1,7 +1,7 @@ /*************************************************************************** * test_GModel.hpp - Test model module * * ----------------------------------------------------------------------- * - * copyright (C) 2012-2016 by Jean-Baptiste Cayrou * + * copyright (C) 2012-2017 by Jean-Baptiste Cayrou * * ----------------------------------------------------------------------- * * * * This program is free software: you can redistribute it and/or modify * @@ -76,6 +76,7 @@ public: void test_spectral_composite(void); void test_spectral_model(void); void test_temp_const(void); + void test_temp_filefct(void); void test_model(void); void test_models(void); void test_model_registry(void); @@ -101,6 +102,7 @@ private: std::string m_map_file; std::string m_cube_file; std::string m_filefct; + std::string m_temp_filefct; std::string m_xml_file; // Model definiton XML files @@ -127,6 +129,7 @@ private: std::string m_xml_model_elliptical_disk; std::string m_xml_model_elliptical_gauss; std::string m_xml_model_spatial_composite; + std::string m_xml_model_point_temp_filefct; // Legacy model definition XML files std::string m_xml_legacy_radial_disk; -- libgit2 0.21.4