From 60616f4d79873c850bfac7af6d8cdf3de9b497bb Mon Sep 17 00:00:00 2001 From: Jürgen Knödlseder Date: Mon, 17 Oct 2016 23:48:17 +0200 Subject: [PATCH] Integrate code --- ChangeLog | 1 + NEWS | 6 ++++++ include/GModelSpectralComposite.hpp | 10 +++++----- pyext/GModelSpectralComposite.i | 9 +++++---- src/model/GModelSpatialComposite.cpp | 2 +- src/model/GModelSpectralComposite.cpp | 222 ++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++------------------------------------------------------------------------------------------ 6 files changed, 150 insertions(+), 100 deletions(-) diff --git a/ChangeLog b/ChangeLog index d2f0f15..7cadfc8 100644 --- a/ChangeLog +++ b/ChangeLog @@ -3,6 +3,7 @@ * Version 1.2.0 released ========================= + Add composite spatial and spectral models (#1706) Add GObservations::save_covmat() method (#1863) Add GFft and GFftWavetable classes (#1731) Add GNdarray class (#1768) diff --git a/NEWS b/NEWS index 3aef732..42c8f4e 100644 --- a/NEWS +++ b/NEWS @@ -15,6 +15,8 @@ The following classes have been added: - GModelSpectralExpInvPlaw - GModelSpectralPlawEnergyFlux - GModelSpectralMultiplicative +- GModelSpectralComposite +- GModelSpatialComposite - GVOTable - GNdarray - GFft @@ -156,6 +158,10 @@ power law. The cut-off parameter is the inverse of the cut-off energy The GModelSpectralMultiplicative class has been added that implements the product of spectral model components as a new spectral model (#1861). +The GModelSpectralComposite and GModelSpatialComposite classes have been +added that implement composite spectral and spatial models, which are sums +of models as a new model component (#1706). + 9. Numerics module ------------------ diff --git a/include/GModelSpectralComposite.hpp b/include/GModelSpectralComposite.hpp index b9f7d57..96a258d 100644 --- a/include/GModelSpectralComposite.hpp +++ b/include/GModelSpectralComposite.hpp @@ -1,7 +1,7 @@ /*************************************************************************** - * GModelSpectralComposite.hpp - Spectral power law model class * + * GModelSpectralComposite.hpp - Spectral composite model class * * ----------------------------------------------------------------------- * - * copyright (C) 2016 by Michael Mayer * + * copyright (C) 2016 by Michael Mayer * * ----------------------------------------------------------------------- * * * * This program is free software: you can redistribute it and/or modify * @@ -82,7 +82,8 @@ public: virtual std::string print(const GChatter& chatter = NORMAL) const; // Other methods - void append(const GModelSpectral& spec, const std::string& name=""); + void append(const GModelSpectral& spec, + const std::string& name=""); int components(void) const; const GModelSpectral* component(const int& index) const; const GModelSpectral* component(const std::string& name) const; @@ -104,7 +105,7 @@ protected: mutable std::vector m_mc_probs; //!< Probailities of individual components mutable GEnergy m_mc_emin; //!< Last minimum energy mutable GEnergy m_mc_emax; //!< Last maximum energy - + mutable std::vector m_mc_values; //!< Parameter values }; @@ -147,5 +148,4 @@ int GModelSpectralComposite::components(void) const return (m_spectral.size()); } - #endif /* GMODELSPECTRALCOMPOSITE_HPP */ diff --git a/pyext/GModelSpectralComposite.i b/pyext/GModelSpectralComposite.i index 3b9d7bc..167c808 100644 --- a/pyext/GModelSpectralComposite.i +++ b/pyext/GModelSpectralComposite.i @@ -1,5 +1,5 @@ /*************************************************************************** - * GModelSpectralComposite.i - Spectral power law model class * + * GModelSpectralComposite.i - Spectral composite model class * * ----------------------------------------------------------------------- * * copyright (C) 2016 Michael Mayer * * ----------------------------------------------------------------------- * @@ -19,7 +19,8 @@ * * ***************************************************************************/ /** - * @file GModelSpecComposite spectral model class interface definition + * @file GModelSpectralComposite.i + * @brief Composite spectral model class interface definition * @author Michael Mayer */ %{ @@ -61,11 +62,11 @@ public: virtual void write(GXmlElement& xml) const; // Other methods - void append(const GModelSpectral& spec, const std::string& name=""); + void append(const GModelSpectral& spec, + const std::string& name=""); int components(void) const; GModelSpectral* component(const int& index) const; GModelSpectral* component(const std::string& name) const; - }; diff --git a/src/model/GModelSpatialComposite.cpp b/src/model/GModelSpatialComposite.cpp index 72b1383..c562790 100644 --- a/src/model/GModelSpatialComposite.cpp +++ b/src/model/GModelSpatialComposite.cpp @@ -662,7 +662,7 @@ void GModelSpatialComposite::free_members(void) for (int i = 0; i < m_components.size(); ++i) { // Delete component i - delete m_components[i]; + if (m_components[i] != NULL) delete m_components[i]; // Signal free pointer m_components[i] = NULL; diff --git a/src/model/GModelSpectralComposite.cpp b/src/model/GModelSpectralComposite.cpp index b0a7336..3401120 100644 --- a/src/model/GModelSpectralComposite.cpp +++ b/src/model/GModelSpectralComposite.cpp @@ -38,15 +38,16 @@ /* __ Constants __________________________________________________________ */ /* __ Globals ____________________________________________________________ */ -const GModelSpectralComposite g_spectral_comp_seed; -const GModelSpectralRegistry g_spectral_comp_registry(&g_spectral_comp_seed); +const GModelSpectralComposite g_spectral_comp_seed; +const GModelSpectralRegistry g_spectral_comp_registry(&g_spectral_comp_seed); /* __ Method name definitions ____________________________________________ */ -#define G_MC "GModelSpectralComposite::mc(GEnergy&, GEnergy&, GTime&, GRan&)" -#define G_WRITE "GModelSpectralComposite::write(GXmlElement&)" -#define G_COMPONENT_INDEX "GModelSpectralComposite::component(const int&)" -#define G_COMPONENT_NAME "GModelSpectralComposite::component(const std::string&)" -#define G_APPEND "GModelSpectralComposite::append(const GModelSpectral&, const std::string&)" +#define G_MC "GModelSpectralComposite::mc(GEnergy&, GEnergy&, GTime&, GRan&)" +#define G_WRITE "GModelSpectralComposite::write(GXmlElement&)" +#define G_COMPONENT_INDEX "GModelSpectralComposite::component(int&)" +#define G_COMPONENT_NAME "GModelSpectralComposite::component(std::string&)" +#define G_APPEND "GModelSpectralComposite::append(GModelSpectral&, "\ + "std::string&)" /* __ Macros _____________________________________________________________ */ @@ -77,13 +78,13 @@ GModelSpectralComposite::GModelSpectralComposite(void) : GModelSpectral() /***********************************************************************//** * @brief XML constructor * - * @param[in] xml XML element containing position information. + * @param[in] xml XML element containing model information. * - * Constructs a power law spectral model by extracting information from an - * XML element. See the read() method for more information about the expected - * structure of the XML element. + * Constructs a composite spectral model by extracting information from an + * XML element. ***************************************************************************/ -GModelSpectralComposite::GModelSpectralComposite(const GXmlElement& xml) : GModelSpectral() +GModelSpectralComposite::GModelSpectralComposite(const GXmlElement& xml) : + GModelSpectral() { // Initialise members init_members(); @@ -99,10 +100,10 @@ GModelSpectralComposite::GModelSpectralComposite(const GXmlElement& xml) : GMode /***********************************************************************//** * @brief Copy constructor * - * @param[in] model Spectral power law model. + * @param[in] model Composite spectral model. ***************************************************************************/ -GModelSpectralComposite::GModelSpectralComposite(const GModelSpectralComposite& model) - : GModelSpectral(model) +GModelSpectralComposite::GModelSpectralComposite(const GModelSpectralComposite& model) : + GModelSpectral(model) { // Initialise members init_members(); @@ -137,8 +138,8 @@ GModelSpectralComposite::~GModelSpectralComposite(void) /***********************************************************************//** * @brief Assignment operator * - * @param[in] model Spectral power law model. - * @return Spectral power law model. + * @param[in] model Composite spectral model. + * @return Composite spectral model. ***************************************************************************/ GModelSpectralComposite& GModelSpectralComposite::operator=(const GModelSpectralComposite& model) { @@ -171,7 +172,7 @@ GModelSpectralComposite& GModelSpectralComposite::operator=(const GModelSpectral ==========================================================================*/ /***********************************************************************//** - * @brief Clear spectral power law model + * @brief Clear composite spectral model ***************************************************************************/ void GModelSpectralComposite::clear(void) { @@ -189,9 +190,9 @@ void GModelSpectralComposite::clear(void) /***********************************************************************//** - * @brief Clone spectral power law model + * @brief Clone composite spectral model * - * @return Pointer to deep copy of spectral power law model. + * @return Pointer to deep copy of composite spectral model. ***************************************************************************/ GModelSpectralComposite* GModelSpectralComposite::clone(void) const { @@ -217,13 +218,12 @@ GModelSpectralComposite* GModelSpectralComposite::clone(void) const * where \f${M_{\rm i}}\f$ is the model evaluation of the i-th model * * If the @p gradients flag is true the method will also compute the - * partial derivatives of each model component with respect to the parameters using - * - * @todo The method expects that energy!=0. Otherwise Inf or NaN may result. + * partial derivatives of each model component with respect to the + * parameters. ***************************************************************************/ double GModelSpectralComposite::eval(const GEnergy& srcEng, - const GTime& srcTime, - const bool& gradients) const + const GTime& srcTime, + const bool& gradients) const { // Initialise result double value = 0.0; @@ -260,11 +260,10 @@ double GModelSpectralComposite::eval(const GEnergy& srcEng, * @param[in] emax Maximum photon energy. * @return Photon flux (ph/cm2/s). * - * Computes the sum of all photon fluxes of each individual component - * + * Computes the sum of the photon fluxes in all components. ***************************************************************************/ double GModelSpectralComposite::flux(const GEnergy& emin, - const GEnergy& emax) const + const GEnergy& emax) const { // Initialise flux double flux = 0.0; @@ -294,11 +293,10 @@ double GModelSpectralComposite::flux(const GEnergy& emin, * @param[in] emax Maximum photon energy. * @return Energy flux (erg/cm2/s). * - * Computes the sum of all energy fluxes of each individual component - * + * Computes the sum of the energy fluxes in all components. ***************************************************************************/ double GModelSpectralComposite::eflux(const GEnergy& emin, - const GEnergy& emax) const + const GEnergy& emax) const { // Initialise eflux double eflux = 0.0; @@ -333,12 +331,13 @@ double GModelSpectralComposite::eflux(const GEnergy& emin, * @exception GException::erange_invalid * Energy range is invalid (emin < emax required). * - * Returns Monte Carlo energy by randomly drawing from a power law. + * Returns Monte Carlo energy by randomly drawing from a composite spectral + * model. ***************************************************************************/ GEnergy GModelSpectralComposite::mc(const GEnergy& emin, - const GEnergy& emax, - const GTime& time, - GRan& ran) const + const GEnergy& emax, + const GTime& time, + GRan& ran) const { // Throw an exception if energy range is invalid if (emin >= emax) { @@ -364,7 +363,6 @@ GEnergy GModelSpectralComposite::mc(const GEnergy& emin, // Loop over spectral components and compute relative probabilites for (int i = 0; i < m_mc_probs.size(); ++i) { - if (u <= m_mc_probs[i]) { index = i; break; @@ -392,7 +390,7 @@ void GModelSpectralComposite::read(const GXmlElement& xml) int n_spectrals = xml.elements("spectrum"); // Loop over spectral elements - for(int i = 0; i < n_spectrals; ++i) { + for (int i = 0; i < n_spectrals; ++i) { // Get spectral XML element const GXmlElement* spec = xml.element("spectrum", i); @@ -495,49 +493,15 @@ void GModelSpectralComposite::write(GXmlElement& xml) const /***********************************************************************//** - * @brief Print powerlaw information - * - * @param[in] chatter Chattiness (defaults to NORMAL). - * @return String containing model information. - ***************************************************************************/ -std::string GModelSpectralComposite::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("=== GModelSpectralComposite ==="); - - // Append information - result.append("\n"+gammalib::parformat("Number of components")); - result.append(gammalib::str(components())); - result.append("\n"+gammalib::parformat("Number of parameters")); - result.append(gammalib::str(size())); - - // Print parameter information - for (int i = 0; i < size(); ++i) { - result.append("\n"+m_pars[i]->print(chatter)); - } - - } // endif: chatter was not silent - - // Return result - return result; -} - - -/***********************************************************************//** * @brief Append spectral component * - * @param[in] spec Spectral model component - * @param[in] name Name of spectral component (can be empty) + * @param[in] spec Spectral model component. + * @param[in] name Name of spectral component. * * Appends a spectral component to the composite model ***************************************************************************/ -void GModelSpectralComposite::append(const GModelSpectral& spec, const std::string& name) +void GModelSpectralComposite::append(const GModelSpectral& spec, + const std::string& name) { // Append model container m_spectral.push_back(spec.clone()); @@ -546,15 +510,15 @@ void GModelSpectralComposite::append(const GModelSpectral& spec, const std::stri int index = m_spectral.size()-1; // Use model index if component name is empty - std::string component_name = !name.empty() ? name : gammalib::str(m_spectral.size()); + std::string component_name = !name.empty() ? name + : gammalib::str(m_spectral.size()); // Check if component name is unique, throw exception if not if (gammalib::contains(m_components, component_name)) { - std::string msg = - "Attempt to append component with name \""+component_name+"\" to composite" - " spectral model container, but a component with the same name exists already.\n" - "Every component in the container needs a unique name. On default the system will" - "increment an integer if no component name is provided."; + std::string msg = "Attempt to append component \""+component_name+"\" " + "to composite spectral model, but a component with " + "the same name exists already. Each component needs " + "a unique name."; throw GException::invalid_value(G_APPEND, msg); } @@ -577,6 +541,9 @@ void GModelSpectralComposite::append(const GModelSpectral& spec, const std::stri m_pars.push_back(par); } // endfor: loop over model parameters + + // Return + return; } @@ -586,18 +553,21 @@ void GModelSpectralComposite::append(const GModelSpectral& spec, const std::stri * @param[in] index Index of spectral component. * @return Spectral model. * + * @exception GException::out_of_range + * Index is out of range. + * * Returns a spectral component to the composite model ***************************************************************************/ const GModelSpectral* GModelSpectralComposite::component(const int& index) const { // Check if index is in validity range if (index >= m_spectral.size() || index < 0) { - throw GException::out_of_range(G_COMPONENT_INDEX, "Component Index", index, m_spectral.size(),""); + throw GException::out_of_range(G_COMPONENT_INDEX, "Component Index", + index, m_spectral.size()); } // Return spectral component return m_spectral[index]; - } @@ -607,6 +577,9 @@ const GModelSpectral* GModelSpectralComposite::component(const int& index) const * @param[in] name Name of spectral component. * @return Spectral model. * + * @exception GException::invalid_argument + * Spectral model not found. + * * Returns a spectral component of the composite model ***************************************************************************/ const GModelSpectral* GModelSpectralComposite::component(const std::string& name) const @@ -622,7 +595,9 @@ const GModelSpectral* GModelSpectralComposite::component(const std::string& name // Check if component name was found if (index == -1) { - throw GException::model_not_found(G_COMPONENT_NAME, name,""); + std::string msg = "Model component \""+name+"\" not found in composite " + "spectral model."; + throw GException::invalid_argument(G_COMPONENT_NAME, msg); } // Return spectral component @@ -630,6 +605,42 @@ const GModelSpectral* GModelSpectralComposite::component(const std::string& name } + +/***********************************************************************//** + * @brief Print composite spectral model information + * + * @param[in] chatter Chattiness. + * @return String containing composite spectral model information. + ***************************************************************************/ +std::string GModelSpectralComposite::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("=== GModelSpectralComposite ==="); + + // Append information + result.append("\n"+gammalib::parformat("Number of components")); + result.append(gammalib::str(components())); + result.append("\n"+gammalib::parformat("Number of parameters")); + result.append(gammalib::str(size())); + + // Print parameter information + 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 = @@ -652,6 +663,7 @@ void GModelSpectralComposite::init_members(void) m_mc_probs.clear(); m_mc_emin.clear(); m_mc_emax.clear(); + m_mc_values.clear(); m_mc_flux = 0.0; // Return @@ -674,6 +686,7 @@ void GModelSpectralComposite::copy_members(const GModelSpectralComposite& model) m_mc_probs = model.m_mc_probs; m_mc_emin = model.m_mc_emin; m_mc_emax = model.m_mc_emax; + m_mc_values = model.m_mc_values; m_mc_flux = model.m_mc_flux; // Copy pointer(s) of spectral component @@ -681,8 +694,10 @@ void GModelSpectralComposite::copy_members(const GModelSpectralComposite& model) // Clear parameters and copy the pointers from the clone m_pars.clear(); - for(int i = 0; i < model.components(); ++i) { - m_spectral.push_back(model.m_spectral[i]->clone()); + for (int i = 0; i < model.components(); ++i) { + + // Clone spectral component + m_spectral.push_back(model.m_spectral[i]->clone()); // Retrieve spectral model GModelSpectral* spec = m_spectral[i]; @@ -711,7 +726,7 @@ void GModelSpectralComposite::copy_members(const GModelSpectralComposite& model) void GModelSpectralComposite::free_members(void) { // Free memory - for(int i = 0; i < m_spectral.size(); ++i) { + for (int i = 0; i < m_spectral.size(); ++i) { // Delete component i if (m_spectral[i] != NULL) delete m_spectral[i]; @@ -734,13 +749,40 @@ void GModelSpectralComposite::free_members(void) * Updates the precomputation cache for Monte Carlo simulations. ***************************************************************************/ void GModelSpectralComposite::update_mc_cache(const GEnergy& emin, - const GEnergy& emax) const + const GEnergy& emax) const { - // Update cache if energy range has changed - if (emin != m_mc_emin || emax != m_mc_emax) { + // Check if one of the parameters has changed. If the dimension of the + // parameter value cache differs from the number of parameters then notify + // a change. This will clear the cache and store the parameter values + // later + bool par_changed = (m_mc_values.size() != size()); + if (par_changed == false) { + for (int i = 0; i < size(); ++i) { + if (m_pars[i]->value() != m_mc_values[i]) { + par_changed = true; + break; + } + } + } + + // Update cache if energy range or parameters have changed + if (par_changed || emin != m_mc_emin || emax != m_mc_emax) { + + // Store energy range m_mc_emin = emin; m_mc_emax = emax; + + // If parameters have changed then store the current parameter + // values for a comparison check for the next method call + if (par_changed) { + m_mc_values.clear(); + for (int i = 0; i < size(); ++i) { + m_mc_values.push_back(m_pars[i]->value()); + } + } + + // Compute MC flux m_mc_flux = flux(m_mc_emin, m_mc_emax); // Initialise sum and probabilites -- libgit2 0.21.4