From 8e14f9ee5f050e07f21173d6c0c607fae05ebfd1 Mon Sep 17 00:00:00 2001 From: Jürgen Knödlseder Date: Tue, 13 Apr 2021 16:37:40 +0200 Subject: [PATCH] Add GModelSpectralBins class (#3606) --- ChangeLog | 3 ++- NEWS | 8 ++++++-- README.md | 4 ++-- doc/source/admin/release_history/2.0.rst | 4 +++- include/GModelSpectralBins.hpp | 178 ++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++ include/GModelSpectralNodes.hpp | 1 - include/GammaLib.hpp | 3 ++- include/Makefile.am | 3 ++- pyext/GModelSpectralBins.i | 107 +++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++ pyext/gammalib/model.i | 3 ++- src/model/GModelSpectralBins.cpp | 1550 ++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++ src/model/Makefile.am | 3 ++- test/data/model_point_bins.xml | 32 ++++++++++++++++++++++++++++++++ test/test_GModel.cpp | 163 ++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++- test/test_GModel.hpp | 4 +++- test/test_GModel.py | 4 ++++ 16 files changed, 2057 insertions(+), 13 deletions(-) create mode 100644 include/GModelSpectralBins.hpp create mode 100644 pyext/GModelSpectralBins.i create mode 100644 src/model/GModelSpectralBins.cpp create mode 100644 test/data/model_point_bins.xml diff --git a/ChangeLog b/ChangeLog index b71c6ac..5ba4984 100644 --- a/ChangeLog +++ b/ChangeLog @@ -1,8 +1,9 @@ -2021-04-02 +2021-04-13 * Version 2.0.0 released ======================== + Add GModelSpectralBins class (#3606) Add carbon footprint estimation to GApplication (#3592) Add phase selection to GCOMSelection and GCOMDri (#3564) Fix incompatibility in GMatrixSparse::mix_column_prepare() method (#3563) diff --git a/NEWS b/NEWS index 6eee9b1..af68a75 100644 --- a/NEWS +++ b/NEWS @@ -1,6 +1,6 @@ New Features and Important Changes in GammaLib 2.0.0 -2 April 2021 +13 April 2021 1. Introduction @@ -19,6 +19,7 @@ The following classes have been added: - GSkyRegionRect - GModelAssociation - GModelAssociations +- GModelSpectralBins - GCTAModelSkyCube - GCTABackground2D @@ -159,8 +160,11 @@ Implemented GMatrixSparse::row() and GMatrixSparse::add_to_row() methods. 8. Model module --------------- +Added GModelSpectralBins class that implements statistically independent spectral +bins as the spectral model component (#3606). + Added methods to GModelSpectralFunc to allow building of a file function -in memort (#3535). +in memory (#3535). Added GModelAssociation class to implement a model association and the GModelAssociations container class to hold a list of model associations. diff --git a/README.md b/README.md index 95fa8d5..7be8545 100644 --- a/README.md +++ b/README.md @@ -1,6 +1,6 @@ GammaLib information ==================== -* Version: 2.0.0.dev (2 April 2021) +* Version: 2.0.0.dev (13 April 2021) [![Build Status](https://cta-jenkins.irap.omp.eu/buildStatus/icon?job=gammalib-integrate-os)](https://cta-jenkins.irap.omp.eu/job/gammalib-integrate-os/) @@ -194,7 +194,7 @@ extensive unit test by typing: If everything works successfully you should see =================== - All 21 tests passed + All 22 tests passed =================== or diff --git a/doc/source/admin/release_history/2.0.rst b/doc/source/admin/release_history/2.0.rst index 9fc0840..f5872f3 100644 --- a/doc/source/admin/release_history/2.0.rst +++ b/doc/source/admin/release_history/2.0.rst @@ -1,6 +1,6 @@ .. _2.0: -GammaLib 2.0 release (13 March 2021) +GammaLib 2.0 release (13 April 2021) ==================================== Overview @@ -73,6 +73,8 @@ Improvements New features ------------ +* [`3606 `_] - + Add ``GModelSpectralBins`` class * [`3527 `_] - Add ``GCTABackground2D`` class * [`3469 `_] - diff --git a/include/GModelSpectralBins.hpp b/include/GModelSpectralBins.hpp new file mode 100644 index 0000000..58fe356 --- /dev/null +++ b/include/GModelSpectralBins.hpp @@ -0,0 +1,178 @@ +/*************************************************************************** + * GModelSpectralBins.hpp - Spectral bins model class * + * ----------------------------------------------------------------------- * + * copyright (C) 2021 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 GModelSpectralBins.hpp + * @brief Spectral bins model class definition + * @author Juergen Knoedlseder + */ + +#ifndef GMODELSPECTRALBINS_HPP +#define GMODELSPECTRALBINS_HPP + +/* __ Includes ___________________________________________________________ */ +#include +#include +#include "GModelSpectral.hpp" +#include "GModelPar.hpp" +#include "GEnergy.hpp" + +/* __ Forward declarations _______________________________________________ */ +class GRan; +class GTime; +class GEbounds; +class GXmlElement; + + +/***********************************************************************//** + * @class GModelSpectralBins + * + * @brief Spectral bins model class + * + * This class implements spectral bins that have a constant intensity within + * their boundaries. A spectral bin is defined by a lower and upper energy + * limit, with the lower limit included and the upper limit excluded from an + * energy bin. + ***************************************************************************/ +class GModelSpectralBins : public GModelSpectral { + +public: + // Constructors and destructors + GModelSpectralBins(void); + GModelSpectralBins(const GModelSpectral& model, + const GEbounds& ebounds, + const double& index = -2.0); + explicit GModelSpectralBins(const GXmlElement& xml); + GModelSpectralBins(const GModelSpectralBins& model); + virtual ~GModelSpectralBins(void); + + // Operators + virtual GModelSpectralBins& operator=(const GModelSpectralBins& model); + + // Implemented pure virtual base class methods + virtual void clear(void); + virtual GModelSpectralBins* clone(void) const; + virtual std::string classname(void) const; + virtual std::string type(void) const; + virtual double eval(const GEnergy& srcEng, + const GTime& srcTime = GTime(), + const bool& gradients = false) const; + virtual double flux(const GEnergy& emin, + const GEnergy& emax) const; + virtual double eflux(const GEnergy& emin, + const GEnergy& emax) const; + virtual GEnergy mc(const GEnergy& emin, + const GEnergy& emax, + const GTime& time, + 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 + int bins(void) const; + void append(const GEnergy& emin, + const GEnergy& emax, + const double& intensity); + void insert(const int& index, + const GEnergy& emin, + const GEnergy& emax, + const double& intensity); + void remove(const int& index); + void reserve(const int& num); + void extend(const GModelSpectralBins& bins); + double index(void) const; + void index(const double& index); + GEnergy emin(const int& index) const; + GEnergy emax(const int& index) const; + void emin(const int& index, const GEnergy& emin); + void emax(const int& index, const GEnergy& emax); + double intensity(const int& index) const; + void intensity(const int& index, const double& intensity); + +protected: + // Protected methods + void init_members(void); + void copy_members(const GModelSpectralBins& model); + void free_members(void); + int bin_index(const GEnergy& energy) const; + void update_pars(void); + void set_cache(void) const; + void mc_update(const GEnergy& emin, const GEnergy& emax) const; + + // Protected members + GModelPar m_index; //!< Spectral index of all bins + std::vector m_emin; //!< Lower energy limits + std::vector m_emax; //!< Upper energy limits + std::vector m_values; //!< Bin values + + // Evaluation cache + mutable std::vector m_epivot; //!< Power-law pivot energies in MeV + + // MC cache + mutable GEnergy m_mc_emin; //!< Minimum energy + mutable GEnergy m_mc_emax; //!< Maximum energy + mutable std::vector m_mc_cum; //!< Cumulative distribution + mutable std::vector m_mc_min; //!< Lower boundary for MC + mutable std::vector m_mc_max; //!< Upper boundary for MC + mutable std::vector m_mc_exp; //!< Exponent for MC +}; + + +/***********************************************************************//** + * @brief Return class name + * + * @return String containing the class name ("GModelSpectralBins"). + ***************************************************************************/ +inline +std::string GModelSpectralBins::classname(void) const +{ + return ("GModelSpectralBins"); +} + + +/***********************************************************************//** + * @brief Return model type + * + * @return "BinFunction". + * + * Returns the type of the spectral bin function model. + ***************************************************************************/ +inline +std::string GModelSpectralBins::type(void) const +{ + return "BinFunction"; +} + + +/***********************************************************************//** + * @brief Return number of bins + * + * @return Number of bins. + * + * Returns the number of bins in the bin function model. + ***************************************************************************/ +inline +int GModelSpectralBins::bins(void) const +{ + return (int)m_values.size(); +} + +#endif /* GMODELSPECTRALBINS_HPP */ diff --git a/include/GModelSpectralNodes.hpp b/include/GModelSpectralNodes.hpp index fe8ec83..9843671 100644 --- a/include/GModelSpectralNodes.hpp +++ b/include/GModelSpectralNodes.hpp @@ -133,7 +133,6 @@ protected: void copy_members(const GModelSpectralNodes& model); void free_members(void); void update_pars(void); - void load_nodes(const std::string& filename); void set_cache(void) const; void set_eval_cache(void) const; void set_flux_cache(void) const; diff --git a/include/GammaLib.hpp b/include/GammaLib.hpp index 31ab9c0..2e75ffa 100644 --- a/include/GammaLib.hpp +++ b/include/GammaLib.hpp @@ -1,7 +1,7 @@ /*************************************************************************** * GammaLib.hpp - Gamma-Ray Astronomy Library Header file * * ----------------------------------------------------------------------- * - * copyright (C) 2006-2020 by Juergen Knoedlseder * + * copyright (C) 2006-2021 by Juergen Knoedlseder * * ----------------------------------------------------------------------- * * * * This program is free software: you can redistribute it and/or modify * @@ -241,6 +241,7 @@ #include "GModelSpectralGauss.hpp" #include "GModelSpectralLogParabola.hpp" #include "GModelSpectralMultiplicative.hpp" +#include "GModelSpectralBins.hpp" #include "GModelSpectralNodes.hpp" #include "GModelSpectralPlaw.hpp" #include "GModelSpectralPlawPhotonFlux.hpp" diff --git a/include/Makefile.am b/include/Makefile.am index 3532f81..f7cbeed 100644 --- a/include/Makefile.am +++ b/include/Makefile.am @@ -1,7 +1,7 @@ ############################################################################# # Makefile.am for include module # # ------------------------------------------------------------------------- # -# Copyright (C) 2016-2020 Juergen Knoedlseder # +# Copyright (C) 2016-2021 Juergen Knoedlseder # # ------------------------------------------------------------------------- # # # # This program is free software: you can redistribute it and/or modify # @@ -200,6 +200,7 @@ pkginclude_HEADERS = GammaLib.hpp \ GModelSpectralGauss.hpp \ GModelSpectralLogParabola.hpp \ GModelSpectralMultiplicative.hpp \ + GModelSpectralBins.hpp \ GModelSpectralNodes.hpp \ GModelSpectralPlaw.hpp \ GModelSpectralPlawPhotonFlux.hpp \ diff --git a/pyext/GModelSpectralBins.i b/pyext/GModelSpectralBins.i new file mode 100644 index 0000000..e4a9e39 --- /dev/null +++ b/pyext/GModelSpectralBins.i @@ -0,0 +1,107 @@ +/*************************************************************************** + * GModelSpectralBins.i - Spectral bins model class * + * ----------------------------------------------------------------------- * + * copyright (C) 2021 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 GModelSpectralBins.i + * @brief Spectral bins model class Python interface definition + * @author Juergen Knoedlseder + */ +%{ +/* Put headers and other declarations here that are needed for compilation */ +#include "GModelSpectralBins.hpp" +%} + + +/***********************************************************************//** + * @class GModelSpectralBins + * + * @brief Spectral bins model class + ***************************************************************************/ +class GModelSpectralBins : public GModelSpectral { +public: + // Constructors and destructors + GModelSpectralBins(void); + GModelSpectralBins(const GModelSpectral& model, + const GEbounds& ebounds, + const double& index = -2.0); + explicit GModelSpectralBins(const GXmlElement& xml); + GModelSpectralBins(const GModelSpectralBins& model); + virtual ~GModelSpectralBins(void); + + // Implemented pure virtual methods + virtual void clear(void); + virtual GModelSpectralBins* clone(void) const; + virtual std::string classname(void) const; + virtual std::string type(void) const; + virtual double eval(const GEnergy& srcEng, + const GTime& srcTime = GTime(), + const bool& gradients = false) const; + virtual double flux(const GEnergy& emin, + const GEnergy& emax) const; + virtual double eflux(const GEnergy& emin, + const GEnergy& emax) const; + virtual GEnergy mc(const GEnergy& emin, + const GEnergy& emax, + const GTime& time, + GRan& ran) const; + virtual void read(const GXmlElement& xml); + virtual void write(GXmlElement& xml) const; + + // Other methods + int bins(void) const; + void append(const GEnergy& emin, + const GEnergy& emax, + const double& intensity); + void insert(const int& index, + const GEnergy& emin, + const GEnergy& emax, + const double& intensity); + void remove(const int& index); + void reserve(const int& num); + void extend(const GModelSpectralBins& bins); + GEnergy emin(const int& index) const; + GEnergy emax(const int& index) const; + void emin(const int& index, const GEnergy& emin); + void emax(const int& index, const GEnergy& emax); + double intensity(const int& index) const; + void intensity(const int& index, const double& intensity); +}; + + +/***********************************************************************//** + * @brief GModelSpectralBins class extension + ***************************************************************************/ +%extend GModelSpectralBins { + GModelSpectralBins copy() { + return (*self); + } +%pythoncode { + def __getstate__(self): + xml = gammalib.GXmlElement() + self.write(xml) + state = xml, + return state + def __setstate__(self, state): + if state[0].elements('bin') == 0: + self.__init__() + else: + self.__init__(state[0]) +} +}; diff --git a/pyext/gammalib/model.i b/pyext/gammalib/model.i index 70831e2..d941f43 100644 --- a/pyext/gammalib/model.i +++ b/pyext/gammalib/model.i @@ -1,7 +1,7 @@ /*************************************************************************** * model.i - Model module * * ----------------------------------------------------------------------- * - * copyright (C) 2011-2020 by Juergen Knoedlseder * + * copyright (C) 2011-2021 by Juergen Knoedlseder * * ----------------------------------------------------------------------- * * * * This program is free software: you can redistribute it and/or modify * @@ -206,6 +206,7 @@ %include "GModelSpectralFunc.i" %include "GModelSpectralGauss.i" %include "GModelSpectralLogParabola.i" +%include "GModelSpectralBins.i" %include "GModelSpectralNodes.i" %include "GModelSpectralPlaw.i" %include "GModelSpectralPlawPhotonFlux.i" diff --git a/src/model/GModelSpectralBins.cpp b/src/model/GModelSpectralBins.cpp new file mode 100644 index 0000000..1d70a4f --- /dev/null +++ b/src/model/GModelSpectralBins.cpp @@ -0,0 +1,1550 @@ +/*************************************************************************** + * GModelSpectralBins.cpp - Spectral bins model class * + * ----------------------------------------------------------------------- * + * copyright (C) 2021 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 GModelSpectralBins.cpp + * @brief Spectral bins model class implementation + * @author Juergen Knoedlseder + */ + +/* __ Includes ___________________________________________________________ */ +#ifdef HAVE_CONFIG_H +#include +#endif +#include +#include "GException.hpp" +#include "GTools.hpp" +#include "GRan.hpp" +#include "GEbounds.hpp" +#include "GModelSpectralBins.hpp" +#include "GModelSpectralRegistry.hpp" + +/* __ Constants __________________________________________________________ */ + +/* __ Globals ____________________________________________________________ */ +const GModelSpectralBins g_spectral_bins_seed; +const GModelSpectralRegistry g_spectral_bins_registry(&g_spectral_bins_seed); + +/* __ Method name definitions ____________________________________________ */ +#define G_FLUX "GModelSpectralBins::flux(GEnergy&, GEnergy&)" +#define G_EFLUX "GModelSpectralBins::eflux(GEnergy&, GEnergy&)" +#define G_MC "GModelSpectralBins::mc(GEnergy&, GEnergy&, GTime&, GRan&)" +#define G_READ "GModelSpectralBins::read(GXmlElement&)" +#define G_WRITE "GModelSpectralBins::write(GXmlElement&)" +#define G_APPEND "GModelSpectralBins::append(GEnergy&, GEnergy&, double&)" +#define G_INSERT "GModelSpectralBins::insert(int&, GEnergy&, GEnergy&, "\ + "double&)" +#define G_REMOVE "GModelSpectralBins::remove(int&)" +#define G_EMIN_GET "GModelSpectralBins::emin(int&)" +#define G_EMAX_GET "GModelSpectralBins::emax(int&)" +#define G_EMIN_SET "GModelSpectralBins::emin(int&, GEnergy&)" +#define G_EMAX_SET "GModelSpectralBins::emax(int&, GEnergy&)" +#define G_INTENSITY_GET "GModelSpectralBins::intensity(int&)" +#define G_INTENSITY_SET "GModelSpectralBins::intensity(int&, double&)" + +/* __ Macros _____________________________________________________________ */ + +/* __ Coding definitions _________________________________________________ */ + +/* __ Debug definitions __________________________________________________ */ + + +/*========================================================================== + = = + = Constructors/destructors = + = = + ==========================================================================*/ + +/***********************************************************************//** + * @brief Void constructor + * + * Constructs an empty spectral bin model. + ***************************************************************************/ +GModelSpectralBins::GModelSpectralBins(void) : GModelSpectral() +{ + // Initialise members + init_members(); + + // Return + return; +} + + +/***********************************************************************//** + * @brief Spectral model constructor + * + * @param[in] model Spectral model. + * @param[in] ebounds Energy boundaries. + * @param[in] index Spectral index. + * + * Constructs a spectral bins model from any spectral model using the + * energy boundaries that are specified by @p ebounds. Within each spectral + * bin the intensity follows a power law with spectral index @p index. + ***************************************************************************/ +GModelSpectralBins::GModelSpectralBins(const GModelSpectral& model, + const GEbounds& ebounds, + const double& index) : + GModelSpectral() +{ + // Initialise members + init_members(); + + // Append bins for all energies + for (int i = 0; i < ebounds.size(); ++i) { + append(ebounds.emin(i), ebounds.emax(i), model.eval(ebounds.elogmean(i))); + } + + // Store index + this->index(index); + + // Return + return; +} + + +/***********************************************************************//** + * @brief XML constructor + * + * @param[in] xml XML element. + * + * Construct spectral bins model by extracting information from an XML + * element. See the read() method for more information about the expected + * structure of the XML element. + ***************************************************************************/ +GModelSpectralBins::GModelSpectralBins(const GXmlElement& xml) : + GModelSpectral() +{ + // Initialise members + init_members(); + + // Read information from XML element + read(xml); + + // Return + return; +} + + +/***********************************************************************//** + * @brief Copy constructor + * + * @param[in] model Spectral bins model. + ***************************************************************************/ +GModelSpectralBins::GModelSpectralBins(const GModelSpectralBins& model) : + GModelSpectral(model) +{ + // Initialise members + init_members(); + + // Copy members + copy_members(model); + + // Return + return; +} + + +/***********************************************************************//** + * @brief Destructor + ***************************************************************************/ +GModelSpectralBins::~GModelSpectralBins(void) +{ + // Free members + free_members(); + + // Return + return; +} + + +/*========================================================================== + = = + = Operators = + = = + ==========================================================================*/ + +/***********************************************************************//** + * @brief Assignment operator + * + * @param[in] model Spectral bins model. + * @return Spectral bins model. + ***************************************************************************/ +GModelSpectralBins& GModelSpectralBins::operator=(const GModelSpectralBins& model) +{ + // Execute only if object is not identical + if (this != &model) { + + // Copy base class members + this->GModelSpectral::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 spectral bins model + ***************************************************************************/ +void GModelSpectralBins::clear(void) +{ + // Free class members (base and derived classes, derived class first) + free_members(); + this->GModelSpectral::free_members(); + + // Initialise members + this->GModelSpectral::init_members(); + init_members(); + + // Return + return; +} + + +/***********************************************************************//** + * @brief Clone spectral bins model +***************************************************************************/ +GModelSpectralBins* GModelSpectralBins::clone(void) const +{ + // Clone spectral bins model + return new GModelSpectralBins(*this); +} + + +/***********************************************************************//** + * @brief Evaluate function + * + * @param[in] srcEng True photon energy. + * @param[in] srcTime True photon arrival time. + * @param[in] gradients Compute gradients? + * @return Model value (ph/cm2/s/MeV). + * + * @todo Document method. + ***************************************************************************/ +double GModelSpectralBins::eval(const GEnergy& srcEng, + const GTime& srcTime, + const bool& gradients) const +{ + // Initialise value + double value = 0.0; + + // Optionally initialise partial derivatives + if (gradients) { + m_index.factor_gradient(0.0); + for (int i = 0; i < m_values.size(); ++i) { + m_values[i].factor_gradient(0.0); + } + } + + // Get bin index + int ibin = bin_index(srcEng); + + // Continue only if index is valid + if (ibin >= 0) { + + // Get power-law index + double index = m_index.value(); + + // Compute and store value + double eng = srcEng.MeV(); + double e_norm = eng / m_epivot[ibin]; + double power = (index != 0.0) ? std::pow(e_norm, index) : 1.0; + + // Compute function value + value = m_values[ibin].value() * power; + + // Compute gradients + if (gradients) { + + // Compute partial derivatives of the parameter values + double g_norm = (m_values[ibin].is_free()) + ? m_values[ibin].scale() * power : 0.0; + double g_index = (m_index.is_free()) + ? value * m_index.scale() * std::log(e_norm) : 0.0; + + // Set gradients + m_values[ibin].factor_gradient(g_norm); + m_index.factor_gradient(g_index); + + } // endif: gradient computation was requested + + // Compile option: Check for NaN/Inf + #if defined(G_NAN_CHECK) + if (gammalib::is_notanumber(value) || gammalib::is_infinite(value)) { + std::cout << "*** ERROR: GModelSpectralBins::eval"; + std::cout << "(srcEng=" << srcEng; + std::cout << ", srcTime=" << srcTime << "):"; + std::cout << " NaN/Inf encountered"; + std::cout << " (value=" << value; + std::cout << ", power=" << power; + std::cout << ")" << std::endl; + } + #endif + + } // endif: bin index was valid + + // Return + return value; +} + + +/***********************************************************************//** + * @brief Returns model photon flux between [emin, emax] (units: ph/cm2/s) + * + * @param[in] emin Minimum photon energy. + * @param[in] emax Maximum photon energy. + * @return Photon flux (ph/cm2/s). + * + * Computes + * + * \f[ + * \int_{\tt emin}^{\tt emax} S_{\rm E}(E | t) dE + * \f] + * + * where + * - [@p emin, @p emax] is an energy interval, and + * - \f$S_{\rm E}(E | t)\f$ is the spectral model (ph/cm2/s/MeV). + ***************************************************************************/ +double GModelSpectralBins::flux(const GEnergy& emin, + const GEnergy& emax) const +{ + // Initialise flux + double flux = 0.0; + + // Compute only if integration range is valid + if (emin < emax) { + + // Get energy range in MeV + double e_min = emin.MeV(); + double e_max = emax.MeV(); + + // Loop over all bins + for (int ibin = 0; ibin < bins(); ++ibin) { + + // Get energy boundaires of bin in MeV + double e_bin_min = m_emin[ibin].value(); + double e_bin_max = m_emax[ibin].value(); + + // Skip bins that do not overlap with energy range + if (e_min >= e_bin_max || e_max < e_bin_min) { + continue; + } + + // Raise lower boundary to minimum energy + if (e_bin_min < e_min) { + e_bin_min = e_min; + } + + // Lower upper boundary to maximum energy + if (e_bin_max > e_max) { + e_bin_max = e_max; + } + + // If energy interval is positive then compute flux + if (e_bin_max > e_bin_min) { + + // Add photon flux + flux += m_values[ibin].value() * + gammalib::plaw_photon_flux(e_bin_min, + e_bin_max, + m_epivot[ibin], + m_index.value()); + + } // endif: energy interval was positive + + } // endfor: looped over all spectral bins + + } // endif: integration range was valid + + // Return + return flux; +} + + +/***********************************************************************//** + * @brief Returns model energy flux between [emin, emax] (units: erg/cm2/s) + * + * @param[in] emin Minimum photon energy. + * @param[in] emax Maximum photon energy. + * @return Energy flux (erg/cm2/s). + * + * Computes + * + * \f[ + * \int_{\tt emin}^{\tt emax} S_{\rm E}(E | t) E \, dE + * \f] + * + * where + * - [@p emin, @p emax] is an energy interval, and + * - \f$S_{\rm E}(E | t)\f$ is the spectral model (ph/cm2/s/MeV). + ***************************************************************************/ +double GModelSpectralBins::eflux(const GEnergy& emin, + const GEnergy& emax) const +{ + // Initialise flux + double eflux = 0.0; + + // Compute only if integration range is valid + if (emin < emax) { + + + // Get energy range in MeV + double e_min = emin.MeV(); + double e_max = emax.MeV(); + + // Loop over all bins + for (int ibin = 0; ibin < bins(); ++ibin) { + + // Get energy boundaires of bin in MeV + double e_bin_min = m_emin[ibin].value(); + double e_bin_max = m_emax[ibin].value(); + + // Skip bins that do not overlap with energy range + if (e_min >= e_bin_max || e_max < e_bin_min) { + continue; + } + + // Raise lower boundary to minimum energy + if (e_bin_min < e_min) { + e_bin_min = e_min; + } + + // Lower upper boundary to maximum energy + if (e_bin_max > e_max) { + e_bin_max = e_max; + } + + // If energy interval is positive then compute flux + if (e_bin_max > e_bin_min) { + + // Add energy flux + eflux += m_values[ibin].value() * + gammalib::plaw_energy_flux(e_bin_min, + e_bin_max, + m_epivot[ibin], + m_index.value()); + + } // endif: energy interval was positive + + } // endfor: looped over all spectral bins + + } // endif: integration range was valid + + // Return flux + return eflux; +} + + +/***********************************************************************//** + * @brief Returns MC energy between [emin, emax] + * + * @param[in] emin Minimum photon energy. + * @param[in] emax Maximum photon energy. + * @param[in] time True photon arrival time. + * @param[in,out] ran Random number generator. + * @return Energy. + * + * @exception GException::erange_invalid + * Energy range is invalid (emin < emax required). + * @exception GException::invalid_return_value + * No valid Monte Carlo cache + * + * Returns Monte Carlo energy by randomly drawing from bin function. + ***************************************************************************/ +GEnergy GModelSpectralBins::mc(const GEnergy& emin, + const GEnergy& emax, + const GTime& time, + GRan& ran) const +{ + // Throw an exception if energy range is invalid + if (emin >= emax) { + throw GException::erange_invalid(G_MC, emin.MeV(), emax.MeV(), + "Minimum energy < maximum energy required."); + } + + // Allocate energy + GEnergy energy; + + // Update cache + mc_update(emin, emax); + + // Determine in which bin we reside. If this fails then thrown an + // exception + 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; + } + } + } + else if (m_mc_cum.size() == 0) { + std::string msg = "No valid bins found for energy interval ["+ + emin.print()+","+emax.print()+"]. Either restrict " + "the energy range to the one covered by the " + "spectral bins or extend the spectral bins " + "in energy."; + throw GException::invalid_return_value(G_MC, msg); + } + + // Get random energy for specific bin + if (m_mc_exp[inx] != 0.0) { + double e_min = m_mc_min[inx]; + double e_max = m_mc_max[inx]; + double u = ran.uniform(); + double eng = (u > 0.0) + ? std::exp(std::log(u * (e_max - e_min) + e_min) / m_mc_exp[inx]) + : 0.0; + energy.MeV(eng); + } + else { + double e_min = m_mc_min[inx]; + double e_max = m_mc_max[inx]; + double u = ran.uniform(); + double eng = std::exp(u * (e_max - e_min) + e_min); + energy.MeV(eng); + } + + // Return energy + return energy; +} + + +/***********************************************************************//** + * @brief Read model from XML element + * + * @param[in] xml XML element. + * + * @exception GException::invalid_value + * Lower energy limit larger than upper energy limit. + * + * Reads the spectral information from an XML element. The format of the XML + * elements is + * + * + * + * + * + * + * + * + * ... + * + * + * + * + * + * + * + * @todo Check that bins are ordered + * @todo Check that energy boundaries are not overlapping + ***************************************************************************/ +void GModelSpectralBins::read(const GXmlElement& xml) +{ + // Free space for bins + m_emin.clear(); + m_emax.clear(); + m_values.clear(); + + // Get index parameter pointer + const GXmlElement* index = gammalib::xml_get_par(G_READ, xml, m_index.name()); + + // Read index parameter + m_index.read(*index); + + // Get number of bins from XML file + int bins = xml.elements("bin"); + + // Loop over all bins + for (int i = 0; i < bins; ++i) { + + // Allocate bin parameters + GModelPar emin; + GModelPar emax; + GModelPar intensity; + + // Get bins + const GXmlElement* bin = xml.element("bin", i); + + // Get parameters + const GXmlElement* lpar = gammalib::xml_get_par(G_READ, *bin, "LowerLimit"); + const GXmlElement* upar = gammalib::xml_get_par(G_READ, *bin, "UpperLimit"); + const GXmlElement* ipar = gammalib::xml_get_par(G_READ, *bin, "Intensity"); + + // Read parameters + emin.read(*lpar); + emax.read(*upar); + intensity.read(*ipar); + + // Throw an exception if either LowerLimit is larger than UpperLimit + if (emin.value() > emax.value()) { + std::string msg = "Lower energy limit "+gammalib::str(emin.value())+ + " MeV is larger than upper energy limit "+ + gammalib::str(emax.value())+" MeV. Please " + "specify lower energy limits that are smaller " + "than the upper energy limits."; + throw GException::invalid_value(G_READ, msg); + } + + // Push bin parameters on list + m_emin.push_back(emin); + m_emax.push_back(emax); + m_values.push_back(intensity); + + } // endfor: looped over bins + + // Update parameter mapping + update_pars(); + + // Set pre-computation cache + set_cache(); + + // Return + return; +} + + +/***********************************************************************//** + * @brief Write model into XML element + * + * @param[in] xml XML element into which model information is written. + * + * @exception GException::invalid_value + * Existing XML element is not of required type + * Invalid number of model parameters or bins found in XML element. + * + * Writes the spectral information into an XML element. The format of the XML + * element is + * + * + * + * + * + * + * + * + * ... + * + * + * + * + * + * + ***************************************************************************/ +void GModelSpectralBins::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 = "Spectral model type \""+xml.attribute("type")+"\" " + "is not of type \""+type()+"\"."; + throw GException::invalid_value(G_WRITE, msg); + } + + // Get index parameter + GXmlElement* index = gammalib::xml_need_par(G_WRITE, xml, m_index.name()); + + // Write index parameter + m_index.write(*index); + + // Determine number of bins + int bins = this->bins(); + + // If XML element has 0 bins then append bins + if (xml.elements("bin") == 0) { + for (int i = 0; i < bins; ++i) { + xml.append(GXmlElement("bin")); + } + } + + // Verify that XML element has the required number of bins and elements. + // Recall that the power-law index also counts as one element. + if (xml.elements() != bins+1) { + std::string msg = "Spectral bins model contains "+ + gammalib::str(xml.elements())+ + " elements but "+gammalib::str(bins+1)+" elements " + "were expected."; + throw GException::invalid_value(G_WRITE, msg); + } + if (xml.elements("bin") != bins) { + std::string msg = "Spectral bins model contains "+ + gammalib::str(xml.elements("bin"))+ + " \"bin\" elements but "+gammalib::str(bins)+ + " \"bin\" elements were expected."; + throw GException::invalid_value(G_WRITE, msg); + } + + // Loop over all bins + for (int i = 0; i < bins; ++i) { + + // Get bins + GXmlElement* bin = xml.element("bin", i); + + // Get copy of parameters so that we can change their names + GModelPar emin = m_emin[i]; + GModelPar emax = m_emax[i]; + GModelPar intensity = m_values[i]; + + // Set names since we appended for internal usage the indices to the + // parameter names, and we want to get rid of them for writing the + // model into the XML element + emin.name("LowerLimit"); + emax.name("UpperLimit"); + intensity.name("Intensity"); + + // Get XML parameters + GXmlElement* lpar = gammalib::xml_need_par(G_WRITE, *bin, emin.name()); + GXmlElement* upar = gammalib::xml_need_par(G_WRITE, *bin, emax.name()); + GXmlElement* ipar = gammalib::xml_need_par(G_WRITE, *bin, intensity.name()); + + // Write parameters + emin.write(*lpar); + emax.write(*upar); + intensity.write(*ipar); + + } // endfor: looped over bins + + // Return + return; +} + + +/***********************************************************************//** + * @brief Append bin + * + * @param[in] emin Lower energy limit. + * @param[in] emax Upper energy limit. + * @param[in] intensity Bin intensity. + * + * @exception GException::invalid_argument + * Lower energy limit is larger than upper energy limit. + * + * Appends one bin to the bin function. By default the lower and upper energy + * limit of the bin is fixed while the intensity of the bin is free. + ***************************************************************************/ +void GModelSpectralBins::append(const GEnergy& emin, + const GEnergy& emax, + const double& intensity) +{ + // Throw an exception if lower limit is larger than upper limit + if (emin > emax) { + std::string msg = "Lower energy limit "+emin.print()+" is larger than " + "upper energy limit "+emax.print()+". Please " + "specify a lower energy limit that is smaller " + "than the upper energy limit."; + throw GException::invalid_value(G_APPEND, msg); + } + + // Allocate bin parameters + GModelPar par_emin; + GModelPar par_emax; + GModelPar par_intensity; + + // Set emin attributes + par_emin.name("LowerLimit"); + par_emin.value(emin.MeV()); + par_emin.unit("MeV"); + par_emin.has_grad(false); + par_emin.fix(); + + // Set emax attributes + par_emax.name("UpperLimit"); + par_emax.value(emax.MeV()); + par_emax.unit("MeV"); + par_emax.has_grad(false); + par_emax.fix(); + + // Set intensity attributes + par_intensity.name("Intensity"); + par_intensity.value(intensity); + par_intensity.unit("ph/cm2/s/MeV"); + par_intensity.has_grad(true); + par_intensity.free(); + + // Append to bins + m_emin.push_back(par_emin); + m_emax.push_back(par_emax); + m_values.push_back(par_intensity); + + // Update parameter mapping + update_pars(); + + // Set pre-computation cache + set_cache(); + + // Return + return; +} + + +/***********************************************************************//** + * @brief Insert bin + * + * @param[in] index Bin index [0,...,bins()-1]. + * @param[in] emin Lower energy limit. + * @param[in] emax Upper energy limit. + * @param[in] intensity Bin intensity. + * + * @exception GException::out_of_range + * Bin index is out of range. + * @exception GException::invalid_argument + * Lower energy limit is larger than upper energy limit. + * + * Inserts a bin into the bin function before the bin with the specified + * @p index. By default the lower and upper energy limit of the bin is fixed + * while the intensity of the bin is free. + ***************************************************************************/ +void GModelSpectralBins::insert(const int& index, + const GEnergy& emin, + const GEnergy& emax, + const double& intensity) +{ + // Raise exception if index is outside boundary + #if defined(G_RANGE_CHECK) + if (index < 0 || index >= bins()) { + throw GException::out_of_range(G_INSERT, "Bin index", index, bins()); + } + #endif + + // Throw an exception if lower limit is larger than upper limit + if (emin > emax) { + std::string msg = "Lower energy limit "+emin.print()+" is larger than " + "upper energy limit "+emax.print()+". Please " + "specify a lower energy limit that is smaller " + "than the upper energy limit."; + throw GException::invalid_value(G_INSERT, msg); + } + + // Allocate bin parameters + GModelPar par_emin; + GModelPar par_emax; + GModelPar par_intensity; + + // Set emin attributes + par_emin.name("LowerLimit"); + par_emin.value(emin.MeV()); + par_emin.unit("MeV"); + par_emin.has_grad(false); + par_emin.fix(); + + // Set emax attributes + par_emax.name("UpperLimit"); + par_emax.value(emax.MeV()); + par_emax.unit("MeV"); + par_emax.has_grad(false); + par_emax.fix(); + + // Set intensity attributes + par_intensity.name("Intensity"); + par_intensity.value(intensity); + par_intensity.unit("ph/cm2/s/MeV"); + par_intensity.has_grad(true); + par_intensity.free(); + + // Insert bin + m_emin.insert(m_emin.begin()+index, par_emin); + m_emax.insert(m_emax.begin()+index, par_emax); + m_values.insert(m_values.begin()+index, par_intensity); + + // Update parameter mapping + update_pars(); + + // Set pre-computation cache + set_cache(); + + // Return + return; +} + + +/***********************************************************************//** + * @brief Remove bin + * + * @param[in] index Bin index [0,...,bins()-1]. + * + * @exception GException::out_of_range + * Bin index is out of range. + * + * Removes bin of specified @p index from the bin function. + ***************************************************************************/ +void GModelSpectralBins::remove(const int& index) +{ + // Raise exception if index is outside boundary + #if defined(G_RANGE_CHECK) + if (index < 0 || index >= bins()) { + throw GException::out_of_range(G_REMOVE, "Bin index", index, bins()); + } + #endif + + // Erase energy and intensity + m_emin.erase(m_emin.begin() + index); + m_emax.erase(m_emax.begin() + index); + m_values.erase(m_values.begin() + index); + + // Update parameter mapping + update_pars(); + + // Set pre-computation cache + set_cache(); + + // Return + return; +} + + +/***********************************************************************//** + * @brief Reserve space for bins + * + * @param[in] num Number of reserved bins. + * + * Reserves space for @p num bins. + ***************************************************************************/ +void GModelSpectralBins::reserve(const int& num) +{ + // Reserve space + m_emin.reserve(num); + m_emax.reserve(num); + m_values.reserve(num); + + // Return + return; +} + + +/***********************************************************************//** + * @brief Append bins from bin function + * + * @param[in] bins Bin function. + * + * Appends all bins from a bin function to current object. + ***************************************************************************/ +void GModelSpectralBins::extend(const GModelSpectralBins& bins) +{ + // Get number of bins in bin function. Note that we extract the size + // first to avoid an endless loop that arises when a container is + // appended to itself. + int num = bins.bins(); + + // Continue only if bin function is not empty + if (num > 0) { + + // Reserve enough space + reserve(this->bins() + num); + + // Loop over all bins and append them to the bin function + for (int i = 0; i < num; ++i) { + m_emin.push_back(bins.m_emin[i]); + m_emax.push_back(bins.m_emax[i]); + m_values.push_back(bins.m_values[i]); + } + + // Update parameter mapping + update_pars(); + + // Set pre-computation cache + set_cache(); + + } // endif: bin function was not empty + + // Return + return; +} + + +/***********************************************************************//** + * @brief Return spectral power law index + * + * @return Spectrel power law index. + * + * Returns the spectral power law index. + ***************************************************************************/ +double GModelSpectralBins::index(void) const +{ + // Return spectral power law index + return (m_index.value()); +} + + +/***********************************************************************//** + * @brief Set spectral power law index + * + * @param[in] index Spectrel power law index. + * + * Set spectral power law index. + ***************************************************************************/ +void GModelSpectralBins::index(const double& index) +{ + // Set m_index + m_index.value(index); + + // Return + return; +} + + +/***********************************************************************//** + * @brief Return lower energy limit of bin + * + * @param[in] index Bin index [0,...,bins()-1]. + * @return Lower energy limit of bin @p index. + * + * @exception GException::out_of_range + * Index is out of range. + * + * Returns the lower energy limit of bin @p index. + ***************************************************************************/ +GEnergy GModelSpectralBins::emin(const int& index) const +{ + // Raise an exception if index is out of range + #if defined(G_RANGE_CHECK) + if (index < 0 || index >= bins()) { + throw GException::out_of_range(G_EMIN_GET, "Bin index", index, bins()); + } + #endif + + // Retrieve energy + GEnergy energy; + energy.MeV(m_emin[index].value()); + + // Return energy + return energy; +} + + +/***********************************************************************//** + * @brief Return upper energy limit of bin + * + * @param[in] index Bin index [0,...,bins()-1]. + * @return Upper energy limit of bin @p index. + * + * @exception GException::out_of_range + * Index is out of range. + * + * Returns the upper energy limit of bin @p index. + ***************************************************************************/ +GEnergy GModelSpectralBins::emax(const int& index) const +{ + // Raise an exception if index is out of range + #if defined(G_RANGE_CHECK) + if (index < 0 || index >= bins()) { + throw GException::out_of_range(G_EMAX_GET, "Bin index", index, bins()); + } + #endif + + // Retrieve energy + GEnergy energy; + energy.MeV(m_emax[index].value()); + + // Return energy + return energy; +} + + +/***********************************************************************//** + * @brief Set lower energy limit of bin + * + * @param[in] index Bin index [0,...,bins()-1]. + * @param[in] emin Lower energy limit of bin @p index. + * + * @exception GException::out_of_range + * Index is out of range. + * + * Sets the lower energy limit of bin @p index. + ***************************************************************************/ +void GModelSpectralBins::emin(const int& index, const GEnergy& emin) +{ + // Raise an exception if index is out of range + #if defined(G_RANGE_CHECK) + if (index < 0 || index >= bins()) { + throw GException::out_of_range(G_EMIN_SET, "Bin index", index, bins()); + } + #endif + + // Set energy + m_emin[index].value(emin.MeV()); + + // Set pre-computation cache + set_cache(); + + // Return + return; +} + + +/***********************************************************************//** + * @brief Set upper energy limit of bin + * + * @param[in] index Bin index [0,...,bins()-1]. + * @param[in] emax Upper energy limit of bin @p index. + * + * @exception GException::out_of_range + * Index is out of range. + * + * Sets the upper energy limit of bin @p index. + ***************************************************************************/ +void GModelSpectralBins::emax(const int& index, const GEnergy& emax) +{ + // Raise an exception if index is out of range + #if defined(G_RANGE_CHECK) + if (index < 0 || index >= bins()) { + throw GException::out_of_range(G_EMAX_SET, "Bin index", index, bins()); + } + #endif + + // Set energy + m_emax[index].value(emax.MeV()); + + // Set pre-computation cache + set_cache(); + + // Return + return; +} + + +/***********************************************************************//** + * @brief Return bin intensity + * + * @param[in] index Bin index [0,...,bins()-1]. + * @return Intensity of bin @p index. + * + * @exception GException::out_of_range + * Index is out of range. + * + * Returns the intensity of bin @p index. + ***************************************************************************/ +double GModelSpectralBins::intensity(const int& index) const +{ + // Raise an exception if index is out of range + #if defined(G_RANGE_CHECK) + if (index < 0 || index >= bins()) { + throw GException::out_of_range(G_INTENSITY_GET, "Bin index", index, bins()); + } + #endif + + // Return intensity + return (m_values[index].value()); +} + + +/***********************************************************************//** + * @brief Set bin intensity + * + * @param[in] index Bin index [0,...,bins()-1]. + * @param[in] intensity Intensity of bin @p index. + * + * @exception GException::out_of_range + * Index is out of range. + * + * Set the intensity of bin @p index. + ***************************************************************************/ +void GModelSpectralBins::intensity(const int& index, const double& intensity) +{ + // Raise an exception if index is out of range + #if defined(G_RANGE_CHECK) + if (index < 0 || index >= bins()) { + throw GException::out_of_range(G_INTENSITY_SET, "Bin index", index, bins()); + } + #endif + + // Set intensity + m_values[index].value(intensity); + + // Return + return; +} + + +/***********************************************************************//** + * @brief Print bin function information + * + * @param[in] chatter Chattiness. + * @return String containing bin function information. + ***************************************************************************/ +std::string GModelSpectralBins::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("=== GModelSpectralBins ==="); + + // Append information + result.append("\n"+gammalib::parformat("Number of bins")); + result.append(gammalib::str(bins())); + 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)); + } + + } // endif: chatter was not silent + + // Return result + return result; +} + + +/*========================================================================== + = = + = Private methods = + = = + ==========================================================================*/ + +/***********************************************************************//** + * @brief Initialise class members + ***************************************************************************/ +void GModelSpectralBins::init_members(void) +{ + // Initialise powerlaw index + m_index.clear(); + m_index.name("Index"); + m_index.scale(1.0); + m_index.value(-2.0); // default: -2.0 + m_index.range(-10.0,+10.0); // range: [-10,+10] + m_index.fix(); + m_index.gradient(0.0); + m_index.has_grad(true); + + // Initialise bin parameter vectors + m_emin.clear(); + m_emax.clear(); + m_values.clear(); + + // Initialise evaluation cache + m_epivot.clear(); + + // Initialise MC cache + m_mc_emin.clear(); + m_mc_emax.clear(); + m_mc_cum.clear(); + m_mc_min.clear(); + m_mc_max.clear(); + m_mc_exp.clear(); + + // Update parameter mapping + update_pars(); + + // Return + return; +} + + +/***********************************************************************//** + * @brief Copy class members + * + * @param[in] model Spectral bin function. + ***************************************************************************/ +void GModelSpectralBins::copy_members(const GModelSpectralBins& model) +{ + // Copy bin energies and values + m_index = model.m_index; + m_emin = model.m_emin; + m_emax = model.m_emax; + m_values = model.m_values; + + // Copy evluation cache + m_epivot = model.m_epivot; + + // Copy MC cache + m_mc_emin = model.m_mc_emin; + m_mc_emax = model.m_mc_emax; + m_mc_cum = model.m_mc_cum; + m_mc_min = model.m_mc_min; + m_mc_max = model.m_mc_max; + m_mc_exp = model.m_mc_exp; + + // Update parameter mapping + update_pars(); + + // Return + return; +} + + +/***********************************************************************//** + * @brief Delete class members + ***************************************************************************/ +void GModelSpectralBins::free_members(void) +{ + // Return + return; +} + + +/***********************************************************************//** + * @brief Return bin index for energy + * + * @param[in] energy Energy. + * @return Bin index (-1 if energy bin was not found). + ***************************************************************************/ +int GModelSpectralBins::bin_index(const GEnergy& energy) const +{ + // Initialise bin index + int bin = -1; + + // Get energy in MeV + double energy_MeV = energy.MeV(); + + // Search bin index + for (int i = 0; i < bins(); ++i) { + if (energy_MeV >= m_emin[i].value() && energy_MeV < m_emax[i].value()) { + bin = i; + break; + } + } + + // Return bin index + return bin; +} + + +/***********************************************************************//** + * @brief Update parameter mapping + * + * Sets the parameter pointers in the m_pars array, enabling iterating over + * all model parameters. This method needs to be called after changing the + * number of bins in the spectral model. The method needs not to be called + * after value update. + ***************************************************************************/ +void GModelSpectralBins::update_pars(void) +{ + // Clear parameter pointer(s) + m_pars.clear(); + + // Push back common spectral index + m_pars.push_back(&m_index); + + // Set parameter pointers for all bins + for (int i = 0; i < bins(); ++i) { + + // Set parameter names + std::string lowerlimit_name = "LowerLimit" + gammalib::str(i); + std::string upperlimit_name = "UpperLimit" + gammalib::str(i); + std::string intensity_name = "Intensity" + gammalib::str(i); + + // Set lower limit attributes + m_emin[i].name(lowerlimit_name); + m_emin[i].unit("MeV"); + m_emin[i].has_grad(false); + + // Set upper limit attributes + m_emax[i].name(upperlimit_name); + m_emax[i].unit("MeV"); + m_emax[i].has_grad(false); + + // Set intensity attributes + m_values[i].name(intensity_name); + m_values[i].unit("ph/cm2/s/MeV"); + m_values[i].has_grad(true); + + // Set pointer + m_pars.push_back(&(m_emin[i])); + m_pars.push_back(&(m_emax[i])); + m_pars.push_back(&(m_values[i])); + + } // endfor: looped over bins + + // Return + return; +} + + +/***********************************************************************//** + * @brief Set cache + * + * Sets the computation cache. The method computes for the time being the + * pivot energies for all bins in MeV. + ***************************************************************************/ +void GModelSpectralBins::set_cache(void) const +{ + // Clear cache + m_epivot.clear(); + + // Compute pivot energies for all bins + for (int i = 0; i < bins(); ++i) { + double epivot = std::sqrt(m_emin[i].value() * m_emax[i].value()); + m_epivot.push_back(epivot); + } + + // Return + return; +} + + +/***********************************************************************//** + * @brief Set MC pre-computation cache + * + * @param[in] emin Minimum energy. + * @param[in] emax Maximum energy. + * + * This method sets up an array of indices and the cumulative distribution + * function needed for MC simulations. + ***************************************************************************/ +void GModelSpectralBins::mc_update(const GEnergy& emin, + const GEnergy& emax) const +{ + // Check if we need to update the cache + if (emin != m_mc_emin || emax != m_mc_emax) { + + // Store new energy interval + m_mc_emin = emin; + m_mc_emax = emax; + + // Initialise cache + m_mc_cum.clear(); + m_mc_min.clear(); + m_mc_max.clear(); + m_mc_exp.clear(); + + // Get energy range in MeV + double e_min = emin.MeV(); + double e_max = emax.MeV(); + + // Continue only if e_max > e_min + if (e_max > e_min) { + + // Loop over all bins + for (int ibin = 0; ibin < bins(); ++ibin) { + + // Get energy boundaires of bin in MeV + double e_bin_min = m_emin[ibin].value(); + double e_bin_max = m_emax[ibin].value(); + + // Skip bins that do not overlap with energy range + if (e_min >= e_bin_max || e_max < e_bin_min) { + continue; + } + + // Raise lower boundary to minimum energy + if (e_bin_min < e_min) { + e_bin_min = e_min; + } + + // Lower upper boundary to maximum energy + if (e_bin_max > e_max) { + e_bin_max = e_max; + } + + // If energy interval is positive then compute flux + if (e_bin_max > e_bin_min) { + + // Get index + double index = m_index.value(); + + // Compute flux + double flux = m_values[ibin].value() * + gammalib::plaw_photon_flux(e_bin_min, + e_bin_max, + m_epivot[ibin], + index); + + // Push results in cache + m_mc_cum.push_back(flux); + m_mc_min.push_back(e_bin_min); + m_mc_max.push_back(e_bin_max); + m_mc_exp.push_back(index); + + } // endif: energy interval was positive + + } // endfor: looped over all spectral bins + + // 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]; + if (norm > 0.0) { + for (int i = 0; i < m_mc_cum.size(); ++i) { + m_mc_cum[i] /= norm; + } + } + + // Set MC values + for (int i = 0; i < m_mc_cum.size(); ++i) { + + // Compute exponent + double exponent = m_mc_exp[i] + 1.0; + + // Exponent dependend computation + if (std::abs(exponent) > 1.0e-11) { + + // If the exponent is too small then use lower energy + // boundary + if (exponent < -50.0) { + m_mc_exp[i] = 0.0; + m_mc_min[i] = std::log(m_mc_min[i]); + m_mc_max[i] = m_mc_min[i]; + } + + // ... otherwise if exponent is too large then use + // upper energy boundary + else if (exponent > +50.0) { + m_mc_exp[i] = 0.0; + m_mc_min[i] = std::log(m_mc_max[i]); + m_mc_max[i] = m_mc_min[i]; + } + + // ... otherwise use transformation formula + else { + m_mc_exp[i] = exponent; + m_mc_min[i] = std::pow(m_mc_min[i], exponent); + m_mc_max[i] = std::pow(m_mc_max[i], exponent); + } + } + else { + m_mc_exp[i] = 0.0; + m_mc_min[i] = std::log(m_mc_min[i]); + m_mc_max[i] = std::log(m_mc_max[i]); + } + + } // endfor: set MC values + + } // endif: e_max > e_min + + } // endif: Update was required + + // Return + return; +} diff --git a/src/model/Makefile.am b/src/model/Makefile.am index 7784065..d2e841e 100644 --- a/src/model/Makefile.am +++ b/src/model/Makefile.am @@ -1,7 +1,7 @@ ############################################################################# # Makefile.am for model module # # ------------------------------------------------------------------------- # -# Copyright (C) 2016-2020 Juergen Knoedlseder # +# Copyright (C) 2016-2021 Juergen Knoedlseder # # ------------------------------------------------------------------------- # # # # This program is free software: you can redistribute it and/or modify # @@ -67,6 +67,7 @@ sources = GModelPar.cpp \ GModelSpectralGauss.cpp \ GModelSpectralLogParabola.cpp \ GModelSpectralMultiplicative.cpp \ + GModelSpectralBins.cpp \ GModelSpectralNodes.cpp \ GModelSpectralPlaw.cpp \ GModelSpectralPlawPhotonFlux.cpp \ diff --git a/test/data/model_point_bins.xml b/test/data/model_point_bins.xml new file mode 100644 index 0000000..5149b59 --- /dev/null +++ b/test/data/model_point_bins.xml @@ -0,0 +1,32 @@ + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + diff --git a/test/test_GModel.cpp b/test/test_GModel.cpp index abde572..2fdd496 100644 --- a/test/test_GModel.cpp +++ b/test/test_GModel.cpp @@ -95,6 +95,7 @@ void TestGModel::set(void) m_xml_model_point_smoothbplaw = datadir + "/model_point_smoothbplaw.xml"; m_xml_model_point_supeplaw = datadir + "/model_point_supeplaw.xml"; m_xml_model_point_logparabola = datadir + "/model_point_logparabola.xml"; + m_xml_model_point_bins = datadir + "/model_point_bins.xml"; m_xml_model_point_nodes = datadir + "/model_point_nodes.xml"; m_xml_model_point_filefct = datadir + "/model_point_filefct.xml"; m_xml_model_point_table = datadir + "/model_point_table.xml"; @@ -201,6 +202,8 @@ void TestGModel::set(void) "Test GModelSpectralExponential"); append(static_cast(&TestGModel::test_nodes), "Test GModelSpectralNodes"); + append(static_cast(&TestGModel::test_bins), + "Test GModelSpectralBins"); append(static_cast(&TestGModel::test_filefct), "Test GModelSpectralFunc"); append(static_cast(&TestGModel::test_table), @@ -2437,6 +2440,164 @@ void TestGModel::test_exponential(void) /***********************************************************************//** + * @brief Test GModelSpectralBins class + ***************************************************************************/ +void TestGModel::test_bins(void) +{ + // Test void constructor + GModelSpectralBins model1; + test_value(model1.type(), "BinFunction", "Check void model type"); + test_value(model1.bins(), 0, "Check number of bins in void model"); + test_value(model1.size(), 1, "Check number of parameters in void model"); + + // Test model constructor + GModelSpectralPlaw plaw(2.0, -2.1, GEnergy(100.0, "MeV")); + GEbounds ebounds(3, GEnergy(1.0, "MeV"), GEnergy(1000.0, "MeV")); + GModelSpectralBins model2(plaw, ebounds, -3.0); + test_value(model2.bins(), 3, "Check number of bins in power-law constructed model"); + test_value(model2.size(), 10, "Check number of parameters in power-law constructed model"); + test_value(model2.index(), -3.0, "Check index in power-law constructed model"); + test_value(model2.emin(0).MeV(), 1.0, "Check emin[0] in power-law constructed model"); + test_value(model2.emax(0).MeV(), 10.0, "Check emax[0] in power-law constructed model"); + test_value(model2.intensity(0), 2825.075089, "Check intensity[0] in power-law constructed model"); + test_value(model2.emin(1).MeV(), 10.0, "Check emin[1] in power-law constructed model"); + test_value(model2.emax(1).MeV(), 100.0, "Check emax[1] in power-law constructed model"); + test_value(model2.intensity(1), 22.440369, "Check intensity[1] in power-law constructed model"); + test_value(model2.emin(2).MeV(), 100.0, "Check emin[2] in power-law constructed model"); + test_value(model2.emax(2).MeV(), 1000.0, "Check emax[2] in power-law constructed model"); + test_value(model2.intensity(2), 0.1782502, "Check intensity[2] in power-law constructed model"); + test_value(model2.eval(GEnergy(1.0, "MeV")), 89336.718430, "Check eval(1 MeV) in power-law constructed model"); + test_value(model2.eval(GEnergy(9.0, "MeV")), 122.546939, "Check eval(9 MeV) in power-law constructed model"); + test_value(model2.eval(GEnergy(76.0, "MeV")), 1.61655029, "Check eval(76 MeV) in power-law constructed model"); + + // Test XML constructor + GXml xml(m_xml_model_point_bins); + GXmlElement* element = xml.element(0)->element(0)->element("spectrum", 0); + GModelSpectralBins model3(*element); + test_value(model3.bins(), 4, "Check number of bins in XML constructed model"); + test_value(model3.size(), 13, "Check number of parameters in XML constructed model"); + test_value(model3.index(), -2.48, "Check index in XML constructed model"); + test_value(model3.emin(0).MeV(), 0.75, "Check emin[0] in XML constructed model"); + test_value(model3.emax(0).MeV(), 1.0, "Check emax[0] in XML constructed model"); + test_value(model3.intensity(0), 1.0e-7, "Check intensity[0] in XML constructed model"); + test_value(model3.emin(1).MeV(), 1.0, "Check emin[1] in XML constructed model"); + test_value(model3.emax(1).MeV(), 3.0, "Check emax[1] in XML constructed model"); + test_value(model3.intensity(1), 0.5e-7, "Check intensity[1] in XML constructed model"); + test_value(model3.emin(2).MeV(), 3.0, "Check emin[2] in XML constructed model"); + test_value(model3.emax(2).MeV(), 10.0, "Check emax[2] in XML constructed model"); + test_value(model3.intensity(2), 0.1e-7, "Check intensity[2] in XML constructed model"); + test_value(model3.emin(3).MeV(), 10.0, "Check emin[3] in XML constructed model"); + test_value(model3.emax(3).MeV(), 30.0, "Check emax[3] in XML constructed model"); + test_value(model3.intensity(3), 0.05e-7, "Check intensity[3] in XML constructed model"); + + // Test copy constructor + GModelSpectralBins model4(model2); + test_value(model4.bins(), 3, "Check number of bins in power-law constructed model"); + test_value(model4.size(), 10, "Check number of parameters in power-law constructed model"); + test_value(model4.index(), -3.0, "Check index in power-law constructed model"); + test_value(model4.emin(0).MeV(), 1.0, "Check emin[0] in power-law constructed model"); + test_value(model4.emax(0).MeV(), 10.0, "Check emax[0] in power-law constructed model"); + test_value(model4.intensity(0), 2825.075089, "Check intensity[0] in power-law constructed model"); + test_value(model4.emin(1).MeV(), 10.0, "Check emin[1] in power-law constructed model"); + test_value(model4.emax(1).MeV(), 100.0, "Check emax[1] in power-law constructed model"); + test_value(model4.intensity(1), 22.440369, "Check intensity[1] in power-law constructed model"); + test_value(model4.emin(2).MeV(), 100.0, "Check emin[2] in power-law constructed model"); + test_value(model4.emax(2).MeV(), 1000.0, "Check emax[2] in power-law constructed model"); + test_value(model4.intensity(2), 0.1782502, "Check intensity[2] in power-law constructed model"); + + // Test XML write and read methods + GXmlElement xml_element; + model4.write(xml_element); + GModelSpectralBins model5; + model5.read(xml_element); + test_value(model5.bins(), 3, "Check number of bins in power-law constructed model"); + + // Test bin function manipulation + GModelSpectralBins model6; + model6.reserve(3); + test_value(model6.bins(), 0, "Check number of bins after first manipulation"); + test_value(model6.size(), 1, "Check number of parameters after first manipulation"); + model6.append(GEnergy(1.0, "MeV"), GEnergy(10.0, "MeV"), 1.0); + test_value(model6.bins(), 1, "Check number of bins after second manipulation"); + test_value(model6.size(), 4, "Check number of parameters after second manipulation"); + test_value(model6.emin(0).MeV(), 1.0, "Check emin[0] after second manipulation"); + test_value(model6.emax(0).MeV(), 10.0, "Check emax[0] after second manipulation"); + test_value(model6.intensity(0), 1.0, "Check intensity[0] after second manipulation"); + model6.append(GEnergy(10.0, "MeV"), GEnergy(100.0, "MeV"), 0.1); + test_value(model6.bins(), 2, "Check number of bins after third manipulation"); + test_value(model6.size(), 7, "Check number of parameters after third manipulation"); + test_value(model6.emin(0).MeV(), 1.0, "Check emin[0] after third manipulation"); + test_value(model6.emax(0).MeV(), 10.0, "Check emax[0] after third manipulation"); + test_value(model6.intensity(0), 1.0, "Check intensity[0] after third manipulation"); + test_value(model6.emin(1).MeV(), 10.0, "Check emin[1] after third manipulation"); + test_value(model6.emax(1).MeV(), 100.0, "Check emax[1] after third manipulation"); + test_value(model6.intensity(1), 0.1, "Check intensity[1] after third manipulation"); + model6.remove(0); + test_value(model6.bins(), 1, "Check number of bins after forth manipulation"); + test_value(model6.size(), 4, "Check number of parameters after forth manipulation"); + test_value(model6.emin(0).MeV(), 10.0, "Check emin[0] after forth manipulation"); + test_value(model6.emax(0).MeV(), 100.0, "Check emax[0] after forth manipulation"); + test_value(model6.intensity(0), 0.1, "Check intensity[0] after forth manipulation"); + model6.insert(0, GEnergy(0.9, "MeV"), GEnergy(9.9, "MeV"), 0.9); + test_value(model6.bins(), 2, "Check number of bins after fifth manipulation"); + test_value(model6.size(), 7, "Check number of parameters after fifth manipulation"); + test_value(model6.emin(0).MeV(), 0.9, "Check emin[0] after fifth manipulation"); + test_value(model6.emax(0).MeV(), 9.9, "Check emax[0] after fifth manipulation"); + test_value(model6.intensity(0), 0.9, "Check intensity[0] after fifth manipulation"); + test_value(model6.emin(1).MeV(), 10.0, "Check emin[1] after fifth manipulation"); + test_value(model6.emax(1).MeV(), 100.0, "Check emax[1] after fifth manipulation"); + test_value(model6.intensity(1), 0.1, "Check intensity[1] after fifth manipulation"); + model6.extend(model6); + test_value(model6.bins(), 4, "Check number of bins after sixth manipulation"); + test_value(model6.size(), 13, "Check number of parameters after sixth manipulation"); + test_value(model6.emin(0).MeV(), 0.9, "Check emin[0] after sixth manipulation"); + test_value(model6.emax(0).MeV(), 9.9, "Check emax[0] after sixth manipulation"); + test_value(model6.intensity(0), 0.9, "Check intensity[0] after sixth manipulation"); + test_value(model6.emin(1).MeV(), 10.0, "Check emin[1] after sixth manipulation"); + test_value(model6.emax(1).MeV(), 100.0, "Check emax[1] after sixth manipulation"); + test_value(model6.intensity(1), 0.1, "Check intensity[1] after sixth manipulation"); + test_value(model6.emin(2).MeV(), 0.9, "Check emin[2] after sixth manipulation"); + test_value(model6.emax(2).MeV(), 9.9, "Check emax[2] after sixth manipulation"); + test_value(model6.intensity(2), 0.9, "Check intensity[2] after sixth manipulation"); + test_value(model6.emin(3).MeV(), 10.0, "Check emin[3] after sixth manipulation"); + test_value(model6.emax(3).MeV(), 100.0, "Check emax[3] after sixth manipulation"); + test_value(model6.intensity(3), 0.1, "Check intensity[3] after sixth manipulation"); + model6.clear(); + test_value(model6.bins(), 0, "Check number of bins after seventh manipulation"); + test_value(model6.size(), 1, "Check number of parameters after seventh manipulation"); + + // Test emin method + model5.emin(0, GEnergy(0.1, "GeV")); + test_value(model5.emin(0).GeV(), 0.1, "Check emin() method"); + + // Test emax method + model5.emax(0, GEnergy(0.2, "GeV")); + test_value(model5.emax(0).GeV(), 0.2, "Check emax() method"); + + // Test intensity method + model5.intensity(0, 2.0e-7); + test_value(model5.intensity(0), 2.0e-7, "Check intensity() method"); + + // Test operator access + const char* strarray[] = {"Index", "LowerLimit0", "UpperLimit0", "Intensity0", + "LowerLimit0", "UpperLimit0", "Intensity1"}; + for (int i = 0; i < 7; ++i) { + std::string keyname(strarray[i]); + model5[keyname].remove_range(); // To allow setting of any value + model5[keyname].value(2.1); + model5[keyname].error(1.9); + model5[keyname].gradient(0.8); + test_value(model5[keyname].value(), 2.1); + test_value(model5[keyname].error(), 1.9); + test_value(model5[keyname].gradient(), 0.8); + } + + // Exit test + return; +} + + +/***********************************************************************//** * @brief Test GModelSpectralNodes class ***************************************************************************/ void TestGModel::test_nodes(void) @@ -2445,7 +2606,7 @@ void TestGModel::test_nodes(void) GModelSpectralNodes model1; test_value(model1.type(), "NodeFunction", "Check void model type"); - // Test node function manimulation + // Test node function manipulation GModelSpectralNodes model2; model2.reserve(3); test_value(model2.size(), 0); diff --git a/test/test_GModel.hpp b/test/test_GModel.hpp index b7a29c6..842b4ae 100644 --- a/test/test_GModel.hpp +++ b/test/test_GModel.hpp @@ -1,7 +1,7 @@ /*************************************************************************** * test_GModel.hpp - Test model module * * ----------------------------------------------------------------------- * - * copyright (C) 2012-2020 by Jean-Baptiste Cayrou * + * copyright (C) 2012-2021 by Jean-Baptiste Cayrou * * ----------------------------------------------------------------------- * * * * This program is free software: you can redistribute it and/or modify * @@ -76,6 +76,7 @@ public: void test_supeplaw(void); void test_logparabola(void); void test_multiplicative(void); + void test_bins(void); void test_nodes(void); void test_filefct(void); void test_spectral_composite(void); @@ -130,6 +131,7 @@ private: std::string m_xml_model_point_smoothbplaw; std::string m_xml_model_point_supeplaw; std::string m_xml_model_point_logparabola; + std::string m_xml_model_point_bins; std::string m_xml_model_point_nodes; std::string m_xml_model_point_filefct; std::string m_xml_model_point_table; diff --git a/test/test_GModel.py b/test/test_GModel.py index c44527c..57dbd53 100755 --- a/test/test_GModel.py +++ b/test/test_GModel.py @@ -143,6 +143,7 @@ class Test(gammalib.GPythonTestSuite): test_support.pickeling(self, gammalib.GModelSpectralGauss()) test_support.pickeling(self, gammalib.GModelSpectralLogParabola()) test_support.pickeling(self, gammalib.GModelSpectralMultiplicative()) + test_support.pickeling(self, gammalib.GModelSpectralBins()) test_support.pickeling(self, gammalib.GModelSpectralNodes()) test_support.pickeling(self, gammalib.GModelSpectralPlaw()) test_support.pickeling(self, gammalib.GModelSpectralPlawEnergyFlux()) @@ -169,6 +170,7 @@ class Test(gammalib.GPythonTestSuite): cnst = gammalib.GModelTemporalConst(3.0) sky = gammalib.GModelSky(ptsrc, plaw, cnst) engs = gammalib.GEnergies(10,pivot,gammalib.GEnergy(10.0,'TeV')) + ebds = gammalib.GEbounds(10,gammalib.GEnergy(1.0,'TeV'),gammalib.GEnergy(10.0,'TeV')) circle = gammalib.GSkyRegionCircle(83.6331, 22.0145, 3.0) spacom = gammalib.GModelSpatialComposite() spacom.append(ptsrc, 'Src1') @@ -181,6 +183,7 @@ class Test(gammalib.GPythonTestSuite): spemul = gammalib.GModelSpectralMultiplicative() spemul.append(plaw, 'Src1') spemul.append(plaw, 'Src2') + spebin = gammalib.GModelSpectralBins(plaw, ebds) spenod = gammalib.GModelSpectralNodes(plaw, engs) spamap = gammalib.GModelSpatialDiffuseMap(cube, 2.0) spamap.mc_cone(circle) @@ -231,6 +234,7 @@ class Test(gammalib.GPythonTestSuite): test_support.pickeling(self, gammalib.GModelSpectralGauss(1.0,pivot,pivot)) test_support.pickeling(self, gammalib.GModelSpectralLogParabola(1.0,-2.0,pivot,0.9)) test_support.pickeling(self, gammalib.GModelSpectralMultiplicative(spemul)) + test_support.pickeling(self, gammalib.GModelSpectralBins(spebin)) test_support.pickeling(self, gammalib.GModelSpectralNodes(spenod)) test_support.pickeling(self, gammalib.GModelSpectralPlaw(plaw)) test_support.pickeling(self, gammalib.GModelSpectralPlawEnergyFlux(1.0,-2.0,pivot,pivot)) -- libgit2 0.21.4