From a6e238ac2fe6708fecc4a1632900134dd240de4e Mon Sep 17 00:00:00 2001 From: knodlseder <> Date: Thu, 17 Feb 2011 22:50:45 +0000 Subject: [PATCH] add model gradient computation method --- include/GObservation.hpp | 46 ++++++++++++++++++++++++++++++++++++++-------- src/obs/GObservation.cpp | 157 +++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++---------------------------- 2 files changed, 167 insertions(+), 36 deletions(-) diff --git a/include/GObservation.hpp b/include/GObservation.hpp index d41ffd4..7ad6e60 100644 --- a/include/GObservation.hpp +++ b/include/GObservation.hpp @@ -96,6 +96,28 @@ protected: void copy_members(const GObservation& obs); void free_members(void); + // Model methods + virtual double model_grad(const GModel& model, const GEvent& event, int ipar) const; + + // Model gradient kernel classes + class model_func : public GFunction { + public: + model_func(const GObservation* parent, + const GModel& model, + const GEvent& event, + int ipar) : + m_parent(parent), + m_model(&model), + m_event(&event), + m_ipar(ipar) { return; } + double eval(double x); + protected: + const GObservation* m_parent; //!< Pointer to parent + const GModel* m_model; //!< Pointer to model + const GEvent* m_event; //!< Pointer to event + int m_ipar; //!< Parameter index + }; + // Npred methods virtual double npred_grad(const GModel& model, int ipar) const; virtual double npred_temp(const GModel& model) const; @@ -104,10 +126,12 @@ protected: // Npred kernel classes class npred_kern_spat : public GIntegrand { public: - npred_kern_spat(const GObservation* parent, const GModel& model, - const GTime& obsTime) : - m_parent(parent), m_model(&model), m_time(&obsTime) - { return; } + npred_kern_spat(const GObservation* parent, + const GModel& model, + const GTime& obsTime) : + m_parent(parent), + m_model(&model), + m_time(&obsTime) { return; } double eval(double x); protected: const GObservation* m_parent; //!< Pointer to parent @@ -116,8 +140,10 @@ protected: }; class npred_kern_spec : public GIntegrand { public: - npred_kern_spec(const GObservation* parent, const GModel& model) : - m_parent(parent), m_model(&model) { return; } + npred_kern_spec(const GObservation* parent, + const GModel& model) : + m_parent(parent), + m_model(&model) { return; } double eval(double x); protected: const GObservation* m_parent; //!< Pointer to parent @@ -127,8 +153,12 @@ protected: // Npred gradient kernel classes class npred_func : public GFunction { public: - npred_func(const GObservation* parent, const GModel& model, int ipar) : - m_parent(parent), m_model(&model), m_ipar(ipar) { ; } + npred_func(const GObservation* parent, + const GModel& model, + int ipar) : + m_parent(parent), + m_model(&model), + m_ipar(ipar) { return; } double eval(double x); protected: const GObservation* m_parent; //!< Pointer to parent diff --git a/src/obs/GObservation.cpp b/src/obs/GObservation.cpp index bb62a8e..bf7c2db 100644 --- a/src/obs/GObservation.cpp +++ b/src/obs/GObservation.cpp @@ -1,5 +1,5 @@ /*************************************************************************** - * GObservation.cpp - Observation abstract base class * + * GObservation.cpp - Abstract observation base class * * ----------------------------------------------------------------------- * * copyright (C) 2008-2011 by Jurgen Knodlseder * * ----------------------------------------------------------------------- * @@ -118,7 +118,7 @@ GObservation::~GObservation(void) * * @param[in] obs Observation. * - * Assign one observation to another. + * Assign observation. ***************************************************************************/ GObservation& GObservation::operator= (const GObservation& obs) { @@ -152,7 +152,7 @@ GObservation& GObservation::operator= (const GObservation& obs) * * @param[in] models Model descriptor. * @param[in] event Observed event. - * @param[out] gradient Pointer to gradient vector (NULL=not computed). + * @param[out] gradient Pointer to gradient vector (optional). * * @exception GException::gradient_par_mismatch * Dimension of gradient vector mismatches number of parameters. @@ -161,21 +161,30 @@ GObservation& GObservation::operator= (const GObservation& obs) * model gives the probability for an event to occur with a given instrument * direction, at a given energy and at a given time. The gradient is the * parameter derivative of this probability. + * + * If NULL is passed for the gradient vector, then gradients will not be + * computed. ***************************************************************************/ double GObservation::model(const GModels& models, const GEvent& event, GVector* gradient) const { - // Verify that gradients vector has the same dimension than the - // model has parameters + // Verify that gradient vector and models have the same dimension #if defined(G_RANGE_CHECK) - if (models.npars() != gradient->size()) - throw GException::gradient_par_mismatch(G_MODEL, gradient->size(), - models.npars()); + if (gradient != NULL) { + if (models.npars() != gradient->size()) + throw GException::gradient_par_mismatch(G_MODEL, + gradient->size(), + models.npars()); + } #endif - // Initialise model value and gradient index - double model = 0.0; - int igrad = 0; + // Initialise + double model = 0.0; // Reset model value + int igrad = 0; // Reset gradient counter + + // If gradient is available then reset gradient vector elements to 0 + if (gradient != NULL) + (*gradient) = 0.0; // Loop over models for (int i = 0; i < models.size(); ++i) { @@ -186,24 +195,18 @@ double GObservation::model(const GModels& models, const GEvent& event, // Compute value and add to model model += models[i].eval_gradients(event, *this); - // Optionally set gradient vector + // Optionally determine model gradients if (gradient != NULL) { - for (int k = 0; k < models[i].size(); ++k, ++igrad) { - double grad = models[i][k].gradient(); - (*gradient)[igrad] = (isinfinite(grad)) ? 0.0 : grad; - } + for (int k = 0; k < models[i].size(); ++k) + (*gradient)[igrad+k] = model_grad(models[i], event, k); } - - } // endif: model was applicable - // ... otherwise set gradient vector to 0 - else if (gradient != NULL) { - for (int k = 0; k < models[i].size(); ++k, ++igrad) { - (*gradient)[igrad] = 0.0; - } - } + } // endif: model component was valid for instrument - } // endfor: looped over models + // Increment parameter counter for gradients + igrad += models[i].size(); + + } // endfor: Looped over models // Return return model; @@ -211,17 +214,31 @@ double GObservation::model(const GModels& models, const GEvent& event, /***********************************************************************//** - * @brief Return total number of predicted counts for all models. + * @brief Return total number (and optionally gradient) of predicted counts + * for all models * * @param[in] models Models. * @param[out] gradient Model parameter gradients (optional). * - * Returns the total number of predicted counts within the analysis region. + * @exception GException::gradient_par_mismatch + * Dimension of gradient vector mismatches number of parameters. * - * @todo Check consistency of gradient verctor + * Returns the total number of predicted counts within the analysis region. + * If NULL is passed for the gradient vector then gradients will not be + * computed. ***************************************************************************/ double GObservation::npred(const GModels& models, GVector* gradient) const { + // Verify that gradient vector and models have the same dimension + #if defined(G_RANGE_CHECK) + if (gradient != NULL) { + if (models.npars() != gradient->size()) + throw GException::gradient_par_mismatch(G_MODEL, + gradient->size(), + models.npars()); + } + #endif + // Initialise double npred = 0.0; // Reset predicted number of counts int igrad = 0; // Reset gradient counter @@ -402,6 +419,48 @@ void GObservation::free_members(void) ==========================================================================*/ /***********************************************************************//** + * @brief Returns parameter gradient of model for a given event + * + * @param[in] model Model. + * @param[in] event Event. + * @param[in] ipar Parameter index for which gradient should be returned. + ***************************************************************************/ +double GObservation::model_grad(const GModel& model, const GEvent& event, + int ipar) const +{ + // Initialise gradient + double grad = 0.0; + + // Compute gradient only if parameter is free + if (model[ipar].isfree()) { + + // If model has a gradient then use it + if (model[ipar].hasgrad()) + grad = model[ipar].gradient(); + + // ... otherwise compute it numerically + else { + + // Setup derivative function + GObservation::model_func function(this, model, event, ipar); + + // Get derivative + GDerivative derivative(&function); + double x = model[ipar].value(); + derivative.eps(1.0e-4); + derivative.silent(true); + grad = derivative.value(x, 0.2*std::abs(x)); + + } // endelse: computed gradient numerically + + } // endif: model parameter was free + + // Return gradient + return grad; +} + + +/***********************************************************************//** * @brief Returns parameter gradient of Npred * * @param[in] model Gamma-ray source model. @@ -625,6 +684,48 @@ double GObservation::npred_kern_spec::eval(double x) /*========================================================================== = = + = Model gradient computation = + = = + ==========================================================================*/ + +/***********************************************************************//** + * @brief Model function evaluation for gradient computation + * + * @param[in] x Function value. + * + * @todo We simply remove any parameter boundaries here for the computation + * to avoid any out of boundary errors. We may have models, however, + * for which out of bound parameters lead to illegal computations, such + * a division by zero or taking the square root of negative values. + * I cannot see any elegant method to catch this at this level. + * Eventually, the higher level method should avoid going in a + * parameter domain that is not defined. + ***************************************************************************/ +double GObservation::model_func::eval(double x) +{ + // Circumvent const correctness + GModel* model = (GModel*)m_model; + + // Save current model parameter + GModelPar current = (*model)[m_ipar]; + + // Set requested model value. Remove any boundaries to avoid limitations. + (*model)[m_ipar].value(x); + (*model)[m_ipar].remove_range(); + + // Compute model value + double value = model->eval_gradients(*m_event, *m_parent); + + // Restore current model parameter + (*model)[m_ipar] = current; + + // Return value + return value; +} + + +/*========================================================================== + = = = Npred gradient integration methods = = = ==========================================================================*/ -- libgit2 0.21.4