Restructure GModel::eval() method for vectorised computations and implement a st…

…ack for model parameters with analytical gradients at the level of the GObservation class.

Restructure GModel::eval() method for vectorised computations and implement a st…
…ack for model parameters with analytical gradients at the level of the GObservation class.
Jürgen Knödlseder
1 parent a133c400
Showing 12 changed files with 351 additions and 193 deletions
include/GModel.hpp
include/GModelData.hpp
include/GModelSky.hpp
include/GObservation.hpp
pyext/GModel.i
pyext/GModelData.i
pyext/GModelSky.i
pyext/GObservation.i
src/model/GModel.cpp
src/model/GModelData.cpp
src/model/GModelSky.cpp
src/obs/GObservation.cpp
@@ -121,16 +121,14 @@ public:
     virtual double      eval(const GEvent& event,
                              const GObservation& obs,
                              const bool& gradients = false) const = 0;
+    virtual GVector     eval(const GObservation& obs,
+                             GMatrixSparse* gradients = NULL) const = 0;
     virtual double      npred(const GEnergy& obsEng, const GTime& obsTime,
                               const GObservation& obs) const = 0;
     virtual void        read(const GXmlElement& xml) = 0;
     virtual void        write(GXmlElement& xml) const = 0;
     virtual std::string print(const GChatter& chatter = NORMAL) const = 0;
  
-    // Virtual methods
-    virtual GVector     eval(const GObservation& obs,
-                             GMatrixSparse* gradients = NULL) const;
-
     // Implemented methods
     int                     size(void) const;
     int                     scales(void) const;
 /***************************************************************************
  *            GModelData.hpp - Abstract virtual data model class           *
  * ----------------------------------------------------------------------- *
- *  copyright (C) 2011-2016 by Juergen Knoedlseder                         *
+ *  copyright (C) 2011-2020 by Juergen Knoedlseder                         *
  * ----------------------------------------------------------------------- *
  *                                                                         *
  *  This program is free software: you can redistribute it and/or modify   *
@@ -40,6 +40,8 @@ class GEvents;
 class GObservation;
 class GRan;
 class GXmlElement;
+class GVector;
+class GMatrixSparse;
  
  
 /***********************************************************************//**
@@ -78,6 +80,10 @@ public:
     virtual void        write(GXmlElement& xml) const = 0;
     virtual std::string print(const GChatter& chatter = NORMAL) const = 0;
  
+    // Implemented pure virtual base class methods
+    virtual GVector     eval(const GObservation& obs,
+                             GMatrixSparse* gradients = NULL) const;
+
 protected:
     // Protected methods
     void init_members(void);
@@ -40,7 +40,6 @@
 #include "GPhoton.hpp"
 #include "GPhotons.hpp"
 #include "GRan.hpp"
-#include "GVector.hpp"
 #include "GEvent.hpp"
 #include "GObservation.hpp"
 #include "GXmlElement.hpp"
@@ -48,6 +47,8 @@
 /* __ Forward declarations _______________________________________________ */
 class GEvent;
 class GObservation;
+class GVector;
+class GMatrixSparse;
  
  
 /***********************************************************************//**
@@ -149,6 +150,8 @@ public:
     virtual double      eval(const GEvent& event,
                              const GObservation& obs,
                              const bool& gradients = false) const;
+    virtual GVector     eval(const GObservation& obs,
+                             GMatrixSparse*      gradients = NULL) const;
     virtual double      npred(const GEnergy& obsEng,
                               const GTime& obsTime,
                               const GObservation& obs) const;
@@ -156,10 +159,6 @@ public:
     virtual void        write(GXmlElement& xml) const;
     virtual std::string print(const GChatter& chatter = NORMAL) const;
  
-    // Reimplemented virtual base class methods
-    virtual GVector     eval(const GObservation& obs,
-                             GMatrixSparse*      gradients = NULL) const;
-
     // Other methods
     GModelSpatial*      spatial(void) const;
     GModelSpectral*     spectral(void) const;
@@ -182,6 +181,7 @@ protected:
     void            free_members(void);
     void            set_pointers(void);
     void            set_type(void);
+    void            signal_analytical_gradients(const GObservation& obs) const;
     GModelSpatial*  xml_spatial(const GXmlElement& spatial) const;
     GModelSpectral* xml_spectral(const GXmlElement& spectral) const;
     GModelTemporal* xml_temporal(const GXmlElement& temporal) const;
@@ -29,17 +29,18 @@
  
 /* __ Includes ___________________________________________________________ */
 #include <string>
+#include <vector>
 #include "GBase.hpp"
 #include "GEvents.hpp"
 #include "GTime.hpp"
 #include "GFunction.hpp"
+#include "GModelPar.hpp"
  
 /* __ Forward declarations _______________________________________________ */
 class GVector;
 class GMatrixSparse;
 class GModel;
 class GModels;
-class GModelPar;
 class GResponse;
 class GXmlElement;
  
@@ -119,12 +120,14 @@ public:
  
     // Implemented methods
     bool               has_events(void) const;
+    bool               has_gradient(const GModelPar& par) const;
     void               name(const std::string& name);
     void               id(const std::string& id);
     void               statistic(const std::string& statistic);
     const std::string& name(void) const;
     const std::string& id(void) const;
     const std::string& statistic(void) const;
+    void               computed_gradient(const GModelPar& par) const;
  
 protected:
     // Protected methods
@@ -217,6 +220,9 @@ protected:
     std::string m_id;          //!< Observation identifier
     std::string m_statistic;   //!< Optimizer statistic
     GEvents*    m_events;      //!< Pointer to event container
+
+    // Stack of pointers to spatial parameters with gradients
+    mutable std::vector<GModelPar*> m_spat_pars_with_gradients;
 };
  
  
@@ -53,6 +53,8 @@ public:
     virtual double      eval(const GEvent& event,
                              const GObservation& obs,
                              const bool& gradients = false) const = 0;
+    virtual GVector     eval(const GObservation& obs,
+                             GMatrixSparse* gradients = NULL) const = 0;
     virtual double      npred(const GEnergy& obsEng, const GTime& obsTime,
                               const GObservation& obs) const = 0;
     virtual void        read(const GXmlElement& xml) = 0;
 /***************************************************************************
  *             GModelData.i - Abstract virtual data model class            *
  * ----------------------------------------------------------------------- *
- *  copyright (C) 2011-2016 by Juergen Knoedlseder                         *
+ *  copyright (C) 2011-2020 by Juergen Knoedlseder                         *
  * ----------------------------------------------------------------------- *
  *                                                                         *
  *  This program is free software: you can redistribute it and/or modify   *
@@ -57,6 +57,10 @@ public:
     virtual GEvents*    mc(const GObservation& obs, GRan& ran) const = 0;
     virtual void        read(const GXmlElement& xml) = 0;
     virtual void        write(GXmlElement& xml) const = 0;
+
+    // Implemented pure virtual base class methods
+    virtual GVector     eval(const GObservation& obs,
+                             GMatrixSparse* gradients = NULL) const;
 };
  
  
 /***************************************************************************
  *                      GModelSky.i - Sky model class                      *
  * ----------------------------------------------------------------------- *
- *  copyright (C) 2011-2018 by Juergen Knoedlseder                         *
+ *  copyright (C) 2011-2020 by Juergen Knoedlseder                         *
  * ----------------------------------------------------------------------- *
  *                                                                         *
  *  This program is free software: you can redistribute it and/or modify   *
@@ -63,6 +63,8 @@ public:
     virtual double      eval(const GEvent& event,
                              const GObservation& obs,
                              const bool& gradients = false) const;
+    virtual GVector     eval(const GObservation& obs,
+                             GMatrixSparse*      gradients = NULL) const;
     virtual double      npred(const GEnergy& obsEng,
                               const GTime& obsTime,
                               const GObservation& obs) const;
 /***************************************************************************
  *             GObservation.i - Abstract observation base class            *
  * ----------------------------------------------------------------------- *
- *  copyright (C) 2008-2019 by Juergen Knoedlseder                         *
+ *  copyright (C) 2008-2020 by Juergen Knoedlseder                         *
  * ----------------------------------------------------------------------- *
  *                                                                         *
  *  This program is free software: you can redistribute it and/or modify   *
@@ -94,12 +94,14 @@ public:
  
     // Implemented methods
     bool               has_events(void) const;
+    bool               has_gradient(const GModelPar& par) const;
     void               name(const std::string& name);
     void               id(const std::string& id);
     void               statistic(const std::string& statistic);
     const std::string& name(void) const;
     const std::string& id(void) const;
     const std::string& statistic(void) const;
+    void               computed_gradient(const GModelPar& par) const;
 };
  
  
@@ -31,13 +31,9 @@
 #include "GTools.hpp"
 #include "GException.hpp"
 #include "GModel.hpp"
-#include "GVector.hpp"
-#include "GMatrixSparse.hpp"
-#include "GObservation.hpp"
  
 /* __ Method name definitions ____________________________________________ */
 #define G_ACCESS                           "GModel::operator[](std::string&)"
-#define G_EVAL                  "GModel::eval(GObservation&, GMatrixSparse*)"
 #define G_AT                                               "GModel::at(int&)"
 #define G_SCALE                              "GModelPar& GModel::scale(int&)"
 #define G_WRITE_SCALES                   "GModel::write_scales(GXmlElement&)"
@@ -225,106 +221,6 @@ const GModelPar&amp; GModel::operator[](const std::string&amp; name) const
  ==========================================================================*/
  
 /***********************************************************************//**
- * @brief Return model values and gradients
- *
- * @param[in] obs Observation.
- * @param[out] gradients Pointer to matrix of gradients.
- * @return Model values.
- *
- * Evaluates the model values and parameter gradients for all events in an
- * observation. Gradients are only returned if the @p gradients pointer is
- * not NULL.
- *
- * The matrix of gradients is a sparse matrix where the number of rows
- * corresponds to the number of model parameters and the number of columns
- * corresponds to the number of events. An exception is thrown if the
- * dimension of the @p gradients matrix is not compatible with the model
- * and the observations.
- ***************************************************************************/
-GVector GModel::eval(const GObservation& obs,
-                     GMatrixSparse*      gradients) const
-{
-    // Get number of model parameters and number of events
-    int npars   = size();
-    int nevents = obs.events()->size();
-
-    // Initialise gradients flag
-    bool grad = (gradients != NULL);
-
-    // Check matrix consistency
-    if (grad) {
-        if (gradients->columns() != npars) {
-            std::string msg = "Number of "+gammalib::str(gradients->columns())+
-                              " columns in gradient matrix differs from number "
-                              "of "+gammalib::str(npars)+" parameters "
-                              "in model. Please provide a compatible gradient "
-                              "matrix.";
-            throw GException::invalid_argument(G_EVAL, msg);
-        }
-        if (gradients->rows() != nevents) {
-            std::string msg = "Number of "+gammalib::str(gradients->rows())+
-                              " rows in gradient matrix differs from number "
-                              "of "+gammalib::str(nevents)+" events in "
-                              "observation. Please provide a compatible "
-                              "gradient matrix.";
-            throw GException::invalid_argument(G_EVAL, msg);
-        }
-    }
-
-    // Initialise values
-    GVector values(nevents);
-
-    // Initialise temporary vectors to hold gradients
-    GVector* tmp_gradients = new GVector[npars];
-    for (int i = 0; i < npars; ++i) {
-        tmp_gradients[i] = GVector(nevents);
-    }
-
-    // Loop over events
-    for (int k = 0; k < nevents; ++k) {
-
-        // Get reference to event
-        const GEvent& event = *((*obs.events())[k]);
-
-        // Evaluate probability
-        values[k] = eval(event, obs, grad);
-
-        // Set gradients
-        if (grad) {
-
-            // Initialise gradient vector
-            //GVector gradient(npars);
-
-            // Extract relevant parameter gradients
-            for (int i = 0; i < npars; ++i) {
-                const GModelPar& par = (*this)[i];
-                if (par.is_free() && par.has_grad()) {
-                    //gradient[i] = par.factor_gradient();
-                    tmp_gradients[i][k] = par.factor_gradient();
-                }
-            }
-
-            // Insert gradient vector into matrix
-            //gradients->add_to_column(k, gradient);
-
-        } // endif: looped over gradients
-
-    } // endfor: looped over events
-
-    // Fill gradients into matrix
-    for (int i = 0; i < npars; ++i) {
-        gradients->column(i, tmp_gradients[i]);
-    }
-
-    // Delete temporal gradients
-    delete [] tmp_gradients;
-
-    // Return values
-    return values;
-}
-
-
-/***********************************************************************//**
  * @brief Returns reference to model parameter by index
  *
  * @param[in] index Parameter index [0,...,size()-1].
 /***************************************************************************
- *           GModelData.cpp  -  Abstract virtual data model class          *
+ *            GModelData.cpp - Abstract virtual data model class           *
  * ----------------------------------------------------------------------- *
- *  copyright (C) 2011 by Jurgen Knodlseder                                *
+ *  copyright (C) 2011-2020 by Juergen Knoedlseder                         *
  * ----------------------------------------------------------------------- *
  *                                                                         *
  *  This program is free software: you can redistribute it and/or modify   *
@@ -21,7 +21,7 @@
 /**
  * @file GModelData.cpp
  * @brief GModelData class implementation.
- * @author J. Knodlseder
+ * @author Juergen Knoedlseder
  */
  
 /* __ Includes ___________________________________________________________ */
@@ -31,8 +31,12 @@
 #include "GTools.hpp"
 #include "GException.hpp"
 #include "GModelData.hpp"
+#include "GVector.hpp"
+#include "GMatrixSparse.hpp"
+#include "GObservation.hpp"
  
 /* __ Method name definitions ____________________________________________ */
+#define G_EVAL                  "GModel::eval(GObservation&, GMatrixSparse*)"
  
 /* __ Macros _____________________________________________________________ */
  
@@ -116,8 +120,9 @@ GModelData::~GModelData(void)
  * @brief Assignment operator
  *
  * @param[in] model Data model.
+ * @return Data model.
  ***************************************************************************/
-GModelData& GModelData::operator= (const GModelData& model)
+GModelData& GModelData::operator=(const GModelData& model)
 {
     // Execute only if object is not identical
     if (this != &model) {
@@ -147,6 +152,120 @@ GModelData&amp; GModelData::operator= (const GModelData&amp; model)
  =                                                                         =
  ==========================================================================*/
  
+/***********************************************************************//**
+ * @brief Return model values and gradients
+ *
+ * @param[in] obs Observation.
+ * @param[out] gradients Pointer to matrix of gradients.
+ * @return Model values.
+ *
+ * Evaluates the model values and parameter gradients for all events in an
+ * observation. Gradients are only returned if the @p gradients pointer is
+ * not NULL.
+ *
+ * The matrix of gradients is a sparse matrix where the number of rows
+ * corresponds to the number of model parameters and the number of columns
+ * corresponds to the number of events. An exception is thrown if the
+ * dimension of the @p gradients matrix is not compatible with the model
+ * and the observations.
+ ***************************************************************************/
+GVector GModelData::eval(const GObservation& obs,
+                         GMatrixSparse*      gradients) const
+{
+    // Get number of model parameters and number of events
+    int npars   = size();
+    int nevents = obs.events()->size();
+
+    // Initialise gradients flag
+    bool grad = ((gradients != NULL) && (npars > 0));
+
+    // Check matrix consistency
+    if (grad) {
+        if (gradients->columns() != npars) {
+            std::string msg = "Number of "+gammalib::str(gradients->columns())+
+                              " columns in gradient matrix differs from number "
+                              "of "+gammalib::str(npars)+" parameters "
+                              "in model. Please provide a compatible gradient "
+                              "matrix.";
+            throw GException::invalid_argument(G_EVAL, msg);
+        }
+        if (gradients->rows() != nevents) {
+            std::string msg = "Number of "+gammalib::str(gradients->rows())+
+                              " rows in gradient matrix differs from number "
+                              "of "+gammalib::str(nevents)+" events in "
+                              "observation. Please provide a compatible "
+                              "gradient matrix.";
+            throw GException::invalid_argument(G_EVAL, msg);
+        }
+    }
+
+    // Initialise values
+    GVector values(nevents);
+
+    // Initialise temporary vectors to hold gradients
+    GVector* tmp_gradients = NULL;
+    if (grad) {
+
+        // Allocate temporary vectors to hold the gradients
+        tmp_gradients = new GVector[npars];
+        for (int i = 0; i < npars; ++i) {
+            tmp_gradients[i] = GVector(nevents);
+        }
+
+        // Signal all parameters that will have analytical gradients. These
+        // are all parameters that are free and for which the model provides
+        // analytical gradients.
+        for (int i = 0; i < npars; ++i) {
+            const GModelPar& par = (*this)[i];
+            if (par.is_free() && par.has_grad()) {
+                obs.computed_gradient(par);
+            }
+        }
+
+    } // endif: gradients were requested
+
+    // Loop over events
+    for (int k = 0; k < nevents; ++k) {
+
+        // Get reference to event
+        const GEvent& event = *((*obs.events())[k]);
+
+        // Evaluate probability
+        values[k] = eval(event, obs, grad);
+
+        // Set gradients
+        if (grad) {
+
+            // Extract relevant parameter gradients
+            for (int i = 0; i < npars; ++i) {
+                const GModelPar& par = (*this)[i];
+                if (par.is_free() && par.has_grad()) {
+                    tmp_gradients[i][k] = par.factor_gradient();
+                }
+            }
+
+        } // endif: looped over gradients
+
+    } // endfor: looped over events
+
+    // Post-process gradients
+    if (grad) {
+
+        // Fill gradients into matrix
+        for (int i = 0; i < npars; ++i) {
+            gradients->column(i, tmp_gradients[i]);
+        }
+
+        // Delete temporal gradients
+        delete [] tmp_gradients;
+
+    } // endif: post-processed gradients
+
+    // Return values
+    return values;
+}
+
+
 /*==========================================================================
  =                                                                         =
  =                             Private methods                             =
@@ -183,10 +302,3 @@ void GModelData::free_members(void)
     // Return
     return;
 }
-
-
-/*==========================================================================
- =                                                                         =
- =                                  Friends                                =
- =                                                                         =
- ==========================================================================*/
@@ -568,6 +568,12 @@ double GModelSky::eval(const GEvent&amp;       event,
     // Evaluate model
     double value = obs.response()->convolve(*this, event, obs, gradients);
  
+    // If gradients were requested then signal all computed analytical
+    // gradients
+    if (gradients) {
+        signal_analytical_gradients(obs);
+    }
+
     // Return
     return value;
 }
@@ -592,6 +598,12 @@ GVector GModelSky::eval(const GObservation&amp; obs,
     // Evaluate model
     GVector values = obs.response()->convolve(*this, obs, gradients);
  
+    // If gradients were requested then signal all computed analytical
+    // gradients
+    if (gradients) {
+        signal_analytical_gradients(obs);
+    }
+
     // Return
     return values;
 }
@@ -1276,6 +1288,63 @@ void GModelSky::set_type(void)
  
  
 /***********************************************************************//**
+ * @brief Signal all parameters that have analytical gradients for a given
+ *        observation
+ *
+ * @param[in] obs Observation.
+ *
+ * This method signals all spectral, temporal and model scaling parameters
+ * for which analytical gradients exist for a given observation. Analytical
+ * gradients will be signalled for all of these parameters that are free
+ * and for which `has_grad()` is `true`.
+ *
+ * Note that spatial analytical gradients are excluded. These will be
+ * signalled by the underlying convolution method as analytical gradients
+ * are not (yet) supported by all convolution methods.
+ ***************************************************************************/
+void GModelSky::signal_analytical_gradients(const GObservation& obs) const
+{
+    // If the model has a spectral component then signal all free parameters
+    // with analytical gradients
+    if (spectral() != NULL) {
+        for (int i = 0; i < spectral()->size(); ++i) {
+            GModelPar& par = (*(spectral()))[i];
+            if (par.is_free() && par.has_grad()) {
+                obs.computed_gradient(par);
+            }
+        }
+    }
+
+    // If the model has a temporal component then signal all free parameters
+    // with analytical gradients
+    if (temporal() != NULL) {
+        for (int i = 0; i < temporal()->size(); ++i) {
+            GModelPar& par = (*(temporal()))[i];
+            if (par.is_free() && par.has_grad()) {
+                obs.computed_gradient(par);
+            }
+        }
+    }
+
+    // If the model has a scales then signal all free parameters with
+    // analytical gradients
+    if (has_scales()) {
+        for (int i = 0; i < scales(); ++i) {
+            const GModelPar& par = scale(i);
+            if (par.name() == obs.instrument()) {
+                if (par.is_free() && par.has_grad()) {
+                    obs.computed_gradient(par);
+                }
+            }
+        }
+    }
+
+    // Return
+    return;
+}
+
+
+/***********************************************************************//**
  * @brief Return pointer to spatial model from XML element
  *
  * @param[in] spatial XML element.
@@ -161,6 +161,72 @@ GObservation&amp; GObservation::operator=(const GObservation&amp; obs)
  ==========================================================================*/
  
 /***********************************************************************//**
+ * @brief Return events
+ *
+ * @exception GException::no_events
+ *            No events allocated for observation.
+ *
+ * Returns pointer to events.
+ ***************************************************************************/
+GEvents* GObservation::events(void)
+{
+    // Throw an exception if the event container is not valid
+    if (m_events == NULL) {
+        std::string msg = "No events allocated for observation.";
+        throw GException::invalid_value(G_EVENTS, msg);
+    }
+
+    // Return pointer to events
+    return (m_events);
+}
+
+
+/***********************************************************************//**
+ * @brief Return events (const version)
+ *
+ * @exception GException::no_events
+ *            No events allocated for observation.
+ *
+ * Returns const pointer to events.
+ ***************************************************************************/
+const GEvents* GObservation::events(void) const
+{
+    // Throw an exception if the event container is not valid
+    if (m_events == NULL) {
+        std::string msg = "No events allocated for observation.";
+        throw GException::invalid_value(G_EVENTS, msg);
+    }
+
+    // Return pointer to events
+    return (m_events);
+}
+
+
+/***********************************************************************//**
+ * @brief Set event container
+ *
+ * @param[in] events Event container.
+ *
+ * Set the event container for this observation by cloning the container
+ * specified in the argument.
+ ***************************************************************************/
+void GObservation::events(const GEvents& events)
+{
+    // Free existing event container only if it differs from current event
+    // container. This prevents unintential deallocation of the argument
+    if ((m_events != NULL) && (m_events != &events)) {
+        delete m_events;
+    }
+
+    // Clone events
+    m_events = events.clone();
+
+    // Return
+    return;
+}
+
+
+/***********************************************************************//**
  * @brief Compute likelihood function
  *
  * @param[in] models Models.
@@ -330,7 +396,7 @@ double GObservation::model(const GModels&amp; models,
  
                             // Set gradient
                             if (par.is_free()) {
-                                if (par.has_grad()) {
+                                if (has_gradient(par)) {
                                     (*gradient)[igrad+ipar] = par.factor_gradient();
                                 }
                                 else {
@@ -356,7 +422,7 @@ double GObservation::model(const GModels&amp; models,
  
                             // Set gradient
                             if (par.is_free()) {
-                                if (par.has_grad()) {
+                                if (has_gradient(par)) {
                                     (*gradient)[igrad+ipar] = par.factor_gradient();
                                 }
                                 else {
@@ -463,11 +529,9 @@ GVector GObservation::model(const GModels&amp; models,
                 else {
  
                     // Allocate gradient matrix
-                    //GMatrixSparse gradients(mptr->size(), nevents);
                     GMatrixSparse gradients(nevents, mptr->size());
  
                     // Initialise sparse matrix stack
-                    //gradients.stack_init(mptr->size()*1000, 1000);
                     gradients.stack_init(nevents*mptr->size(), mptr->size());
  
                     // Evaluate model and add to values
@@ -495,14 +559,11 @@ GVector GObservation::model(const GModels&amp; models,
  
                             // Determine gradient
                             GVector grad(nevents);
-                            if (par.is_free()) {
-                                if (par.has_grad()) {
-                                    //grad = gradients.row(ipar);
-                                    grad = gradients.column(ipar);
-                                }
-                                else {
-                                    grad = model_grad(*mptr, par);
-                                }
+                            if (has_gradient(par)) {
+                                grad = gradients.column(ipar);
+                            }
+                            else {
+                                grad = model_grad(*mptr, par);
                             }
  
                             // Set gradient
@@ -523,10 +584,8 @@ GVector GObservation::model(const GModels&amp; models,
  
                             // Determine gradient
                             GVector grad(nevents);
-                            if (par.is_free()) {
-                                if (par.has_grad()) {
-                                    //grad = gradients.row(ipar);
-                                    grad = gradients.column(ipar);
+                            if (has_gradient(par)) {
+                                grad = gradients.column(ipar);
 // Debugging
 /*
 if (ipar < 3) {
@@ -542,10 +601,9 @@ if (ipar &lt; 3) {
     }
 }
 */
-                                }
-                                else {
-                                    grad = model_grad(*mptr, par);
-                                }
+                            }
+                            else {
+                                grad = model_grad(*mptr, par);
                             }
  
                             // Set gradient
@@ -1110,68 +1168,65 @@ void GObservation::remove_response_cache(const std::string&amp; name)
  
  
 /***********************************************************************//**
- * @brief Set event container
+ * @brief Check whether a model parameter has an analytical gradient
  *
- * @param[in] events Event container.
+ * @param[in] par Model parameter.
+ * @returns True if model parameter is free and has an analytical gradient.
  *
- * Set the event container for this observation by cloning the container
- * specified in the argument.
+ * Checks whether a model parameter is free and has an analytical gradient.
  ***************************************************************************/
-void GObservation::events(const GEvents& events)
+bool GObservation::has_gradient(const GModelPar& par) const
 {
-    // Free existing event container only if it differs from current event
-    // container. This prevents unintential deallocation of the argument
-    if ((m_events != NULL) && (m_events != &events)) {
-        delete m_events;
-    }
+    // Initialise flag
+    bool has_gradient = false;
+    
+    // Only consider free parameters
+    if (par.is_free()) {
  
-    // Clone events
-    m_events = events.clone();
+        // Search for parameter address in stack and set flag to true if
+        // parameter address was found
+        for (int i = 0; i < m_spat_pars_with_gradients.size(); ++i) {
+            if (m_spat_pars_with_gradients[i] == &par) {
+                has_gradient = true;
+                break;
+            }
+        }
  
-    // Return
-    return;
+    } // endif: parameter was free
+
+    // Return flag
+    return has_gradient;
 }
  
  
 /***********************************************************************//**
- * @brief Return events
+ * @brief Signals that an analytical gradient was computed for a model
+ *        parameter
  *
- * @exception GException::no_events
- *            No events allocated for observation.
+ * @param[in] par Model parameter.
  *
- * Returns pointer to events.
+ * Signals that an analytical gradient was computed for a model parameter.
  ***************************************************************************/
-GEvents* GObservation::events(void)
+void GObservation::computed_gradient(const GModelPar& par) const
 {
-    // Throw an exception if the event container is not valid
-    if (m_events == NULL) {
-        std::string msg = "No events allocated for observation.";
-        throw GException::invalid_value(G_EVENTS, msg);
+    // Initialise flag
+    bool in_stack = false;
+
+    // Check if parameter is already in stack
+    for (int i = 0; i < m_spat_pars_with_gradients.size(); ++i) {
+        if (m_spat_pars_with_gradients[i] == &par) {
+            in_stack = true;
+            break;
+        }
     }
  
-    // Return pointer to events
-    return (m_events);
-}
-
-
-/***********************************************************************//**
- * @brief Return events (const version)
- *
- * @exception GException::no_events
- *            No events allocated for observation.
- *
- * Returns const pointer to events.
- ***************************************************************************/
-const GEvents* GObservation::events(void) const
-{
-    // Throw an exception if the event container is not valid
-    if (m_events == NULL) {
-        std::string msg = "No events allocated for observation.";
-        throw GException::invalid_value(G_EVENTS, msg);
+    // If parameter is not yet in stack the push it on the stack
+    if (!in_stack) {
+        m_spat_pars_with_gradients.push_back(const_cast<GModelPar*>(&par));
     }
  
-    // Return pointer to events
-    return (m_events);
+    // Return
+    return;
 }
  
  
@@ -1192,6 +1247,9 @@ void GObservation::init_members(void)
     m_statistic = "cstat";
     m_events    = NULL;
  
+    // Initialise stack of pointers to spatial parameters with gradients
+    m_spat_pars_with_gradients.clear();
+
     // Return
     return;
 }
@@ -1211,6 +1269,9 @@ void GObservation::copy_members(const GObservation&amp; obs)
     m_id        = obs.m_id;
     m_statistic = obs.m_statistic;
  
+    // Copy stack of pointers to spatial parameters with gradients
+    m_spat_pars_with_gradients = obs.m_spat_pars_with_gradients;
+
     // Clone members
     m_events = (obs.m_events != NULL) ? obs.m_events->clone() : NULL;