Integrate code

Jürgen Knödlseder
1 parent 10f11c99
Showing 6 changed files with 150 additions and 100 deletions
ChangeLog
NEWS
include/GModelSpectralComposite.hpp
pyext/GModelSpectralComposite.i
src/model/GModelSpatialComposite.cpp
src/model/GModelSpectralComposite.cpp
@@ -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)
@@ -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
 ------------------
 /***************************************************************************
- *         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<double>  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<double>  m_mc_values;   //!< Parameter values
 };
  
  
@@ -147,5 +148,4 @@ int GModelSpectralComposite::components(void) const
     return (m_spectral.size());
 }
  
-
 #endif /* GMODELSPECTRALCOMPOSITE_HPP */
 /***************************************************************************
- *          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;
-
 };
  
  
@@ -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;
@@ -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&amp; 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&amp; 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&amp; 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&amp; 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&amp; 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&amp; 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&amp; 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&amp; 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&amp; 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&amp; 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&amp; 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&amp; 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&amp; 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&amp; 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&amp; 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&amp; 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&amp; 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