From d9f99bde0291386964b345817bc5eb16ff600fed Mon Sep 17 00:00:00 2001 From: Jürgen Knödlseder Date: Wed, 22 May 2019 21:07:28 +0200 Subject: [PATCH] Allow fitting of instrument scale factors (#2892) --- ChangeLog | 3 ++- NEWS | 7 ++++++- doc/source/admin/release_history/1.6.rst | 2 ++ include/GModel.hpp | 19 ++++++++++++++++++- inst/cta/src/GCTAOnOffObservation.cpp | 47 ++++++++++++++++++++++++++++++++++++----------- pyext/GModel.i | 8 +++++--- src/model/GModel.cpp | 69 +++++++++++++++++++++++++++++++++++++++++++++++++++------------------ src/model/GModelSky.cpp | 15 ++++++++++++++- src/obs/GObservation.cpp | 2 +- src/obs/GResponse.cpp | 35 +++++++++++++++++++++++++++-------- test/test_GModel.cpp | 7 +++++-- 11 files changed, 167 insertions(+), 47 deletions(-) diff --git a/ChangeLog b/ChangeLog index 6ff9e21..2d57d80 100644 --- a/ChangeLog +++ b/ChangeLog @@ -1,8 +1,9 @@ -2019-04-26 +2019-05-22 * Version 1.6.0 released ======================== + Allow fitting of instrument scale factors (#2892) Add GCTAModelSpatialLookup class to handle CTA background lookups (#2878) Implement energy-dependent Gaussian CTA spatial background model (#2868) Take instrument scale factors in On/Off analysis into account (#2867) diff --git a/NEWS b/NEWS index f9c300e..8d0b2e5 100644 --- a/NEWS +++ b/NEWS @@ -1,6 +1,6 @@ New Features and Important Changes in GammaLib 1.6.0 -26 April 2019 +22 May 2019 1. Introduction @@ -73,6 +73,8 @@ The following methods have been added: - GFitsTable::remove(int&) - GFitsTable::remove(std::string&) - GObservation::has_events() +- GModel::scales() +- GModel::scale(int&) - GModelSpatialPointSource::type(std::string&) - GModelSpectralPlaw::type(std::string&) - GModelSpatial::has_free_pars() @@ -217,6 +219,9 @@ None 8. Model module --------------- +Add GModel::scales() method to return number of scales for this model and +add GModel::scale() parameter access methods (#2892). + Add GModelSpatial::has_free_pars() method that checks whether a spatial model component has free model parameters or not (#2715). diff --git a/doc/source/admin/release_history/1.6.rst b/doc/source/admin/release_history/1.6.rst index 53dbec8..bdfb8d9 100644 --- a/doc/source/admin/release_history/1.6.rst +++ b/doc/source/admin/release_history/1.6.rst @@ -55,6 +55,8 @@ Bug fixes Improvements ------------ +* [`2892 `_] - + Allow fitting of instrument scale factors * [`2878 `_] - Add ``GCTAModelSpatialLookup`` class to handle CTA background lookups * [`2868 `_] - diff --git a/include/GModel.hpp b/include/GModel.hpp index ee64ca3..29246df 100644 --- a/include/GModel.hpp +++ b/include/GModel.hpp @@ -1,7 +1,7 @@ /*************************************************************************** * GModel.hpp - Abstract virtual model base class * * ----------------------------------------------------------------------- * - * copyright (C) 2009-2016 by Juergen Knoedlseder * + * copyright (C) 2009-2019 by Juergen Knoedlseder * * ----------------------------------------------------------------------- * * * * This program is free software: you can redistribute it and/or modify * @@ -127,6 +127,7 @@ public: // Implemented methods int size(void) const; + int scales(void) const; GModelPar& at(const int& index); const GModelPar& at(const int& index) const; bool has_par(const std::string& name) const; @@ -140,6 +141,8 @@ public: const bool& has_ts(void) const; std::string instruments(void) const; void instruments(const std::string& instruments); + GModelPar& scale(const int& index); + const GModelPar& scale(const int& index) const; GModelPar scale(const std::string& instrument) const; void scale(const GModelPar& par); std::string ids(void) const; @@ -218,6 +221,20 @@ int GModel::size(void) const /***********************************************************************//** + * @brief Return number of scale parameters in model + * + * @return Number of scale parameters in model. + * + * Returns the number of scale parameters in the model. + ***************************************************************************/ +inline +int GModel::scales(void) const +{ + return (int)m_scales.size(); +} + + +/***********************************************************************//** * @brief Return parameter name * * @return Parameter name. diff --git a/inst/cta/src/GCTAOnOffObservation.cpp b/inst/cta/src/GCTAOnOffObservation.cpp index 0e5c778..17b624d 100644 --- a/inst/cta/src/GCTAOnOffObservation.cpp +++ b/inst/cta/src/GCTAOnOffObservation.cpp @@ -2329,6 +2329,16 @@ double GCTAOnOffObservation::N_gamma(const GModels& models, double rmf_sum = 0.0; #endif + // Get indices of spectral gradients. We need this to update + // only the spectral parameter gradients later. + std::vector inx_spec; + for (int k = 0; k < mptr->size(); ++k) { + const GModelPar& par = (*mptr)[k]; + if (spectral->has_par(par.name())) { + inx_spec.push_back(k); + } + } + // Set instrument scale factor double scale = (sky->has_scales()) ? sky->scale(instrument()).value() : 1.0; @@ -2358,16 +2368,20 @@ double GCTAOnOffObservation::N_gamma(const GModels& models, GEnergy etruemean = m_arf.ebounds().elogmean(itrue); double etruewidth = m_arf.ebounds().ewidth(itrue).MeV(); - // Compute normalisation factors - double exposure = m_on_spec.exposure(); - double norm_flux = arf * exposure * rmf * scale; - double norm_grad = norm_flux * etruewidth; + // Compute normalisation factors. The gradient normalisation + // factor uses the energy bin width, since the eval() method + // returns a differential flux + double exposure = m_on_spec.exposure(); + double norm_scale = arf * exposure * rmf; + double norm_flux = norm_scale * scale; + double norm_grad = norm_flux * etruewidth; // Determine number of gamma-ray events in model by // computing the flux over the true energy bin in // ph/cm2/s and multiplying this by effective area (cm2), - // livetime (s) and redistribution probability - value += spectral->flux(etruemin, etruemax) * norm_flux; + // livetime (s), redistribution probability and scale + double flux = spectral->flux(etruemin, etruemax); + value += flux * norm_flux; // Determine the model gradients at the current true // energy. The eval() method needs a time in case that the @@ -2375,11 +2389,22 @@ double GCTAOnOffObservation::N_gamma(const GModels& models, // dummy time here. spectral->eval(etruemean, GTime(), true); - // Compute the parameter gradients for all model parameters - for (int k = 0; k < mptr->size(); ++k) { - const GModelPar& par = (*mptr)[k]; - if (par.is_free() && ipar+k < npars) { - (*grad)[ipar+k] += par.factor_gradient() * norm_grad; + // Compute the parameter gradients for spectral parameters + for (int k = 0; k < inx_spec.size(); ++k) { + const GModelPar& par = (*mptr)[inx_spec[k]]; + if (par.is_free() && ipar+k < npars) { + (*grad)[ipar+inx_spec[k]] += par.factor_gradient() * norm_grad; + } + } + + // Optionally compute scaling parameter gradient + if (sky->has_scales()) { + for (int k = 0; k < mptr->size(); ++k) { + const GModelPar& par = (*mptr)[k]; + if ((par.name() == instrument()) && par.is_free() && + (ipar+k < npars)) { + (*grad)[ipar+k] += par.scale() * flux * norm_scale; + } } } diff --git a/pyext/GModel.i b/pyext/GModel.i index 13dd1b2..9c4d3ac 100644 --- a/pyext/GModel.i +++ b/pyext/GModel.i @@ -1,7 +1,7 @@ /*************************************************************************** * GModel.i - Abstract virtual model base class * * ----------------------------------------------------------------------- * - * copyright (C) 2009-2018 by Juergen Knoedlseder * + * copyright (C) 2009-2019 by Juergen Knoedlseder * * ----------------------------------------------------------------------- * * * * This program is free software: you can redistribute it and/or modify * @@ -60,17 +60,19 @@ public: // Implemented methods int size(void) const; + int scales(void) const; bool has_par(const std::string& name) const; bool has_scales(void) const; const std::string& name(void) const; void name(const std::string& name); - std::string instruments(void) const; - void instruments(const std::string& instruments); const double& ts(void) const; void ts(const double& ts); const bool& tscalc(void) const; void tscalc(const bool& tscalc); const bool& has_ts(void) const; + std::string instruments(void) const; + void instruments(const std::string& instruments); + GModelPar& scale(const int& index); GModelPar scale(const std::string& instrument) const; void scale(const GModelPar& par); std::string ids(void) const; diff --git a/src/model/GModel.cpp b/src/model/GModel.cpp index 9257290..d55aebf 100644 --- a/src/model/GModel.cpp +++ b/src/model/GModel.cpp @@ -1,7 +1,7 @@ /*************************************************************************** * GModel.cpp - Abstract virtual model base class * * ----------------------------------------------------------------------- * - * copyright (C) 2009-2015 by Juergen Knoedlseder * + * copyright (C) 2009-2019 by Juergen Knoedlseder * * ----------------------------------------------------------------------- * * * * This program is free software: you can redistribute it and/or modify * @@ -35,6 +35,7 @@ /* __ Method name definitions ____________________________________________ */ #define G_ACCESS "GModel::operator[](std::string&)" #define G_AT "GModel::at(int&)" +#define G_SCALE "GModelPar& GModel::scale(int&)" #define G_WRITE_SCALES "GModel::write_scales(GXmlElement&)" /* __ Macros _____________________________________________________________ */ @@ -348,6 +349,54 @@ void GModel::instruments(const std::string& instruments) /***********************************************************************//** + * @brief Returns reference to scale parameter by index + * + * @param[in] index Parameter index [0,...,size()-1]. + * @return Reference to scale parameter. + * + * @exception GException::out_of_range + * Scale parameter index is out of range. + * + * Returns a reference to the scale parameter of the specified @p index. + * Throws an exception if @p index is not valid. + ***************************************************************************/ +GModelPar& GModel::scale(const int& index) +{ + // Raise exception if index is out of range + if (index < 0 || index >= scales()) { + throw GException::out_of_range(G_SCALE, index, scales()); + } + + // Return reference + return (m_scales[index]); +} + + +/***********************************************************************//** + * @brief Returns reference to scale parameter by index (const version) + * + * @param[in] index Parameter index [0,...,size()-1]. + * @return Reference to scale parameter. + * + * @exception GException::out_of_range + * Scale parameter index is out of range. + * + * Returns a reference to the scale parameter of the specified @p index. + * Throws an exception if @p index is not valid. + ***************************************************************************/ +const GModelPar& GModel::scale(const int& index) const +{ + // Raise exception if index is out of range + if (index < 0 || index >= scales()) { + throw GException::out_of_range(G_SCALE, index, scales()); + } + + // Return reference + return (m_scales[index]); +} + + +/***********************************************************************//** * @brief Returns model scale factor for a given instrument * * @param[in] instrument Instrument. @@ -715,23 +764,6 @@ std::string GModel::print_attributes(void) const result.append("Computation requested"); } - // Append instrument scale factors - result.append("\n"+gammalib::parformat("Instrument scale factors")); - if (!m_scales.empty()) { - for (int i = 0; i < m_scales.size(); ++i) { - if (i > 0) { - result.append(", "); - } - result.append(m_scales[i].name()); - result.append("="); - result.append(gammalib::str(m_scales[i].value())); - } - result.append(", others unity"); - } - else { - result.append("unity"); - } - // Append observation identifiers result.append("\n"+gammalib::parformat("Observation identifiers")); if (!m_ids.empty()) { @@ -788,6 +820,7 @@ void GModel::read_scales(const GXmlElement& xml) GModelPar scale; scale.read(*par); scale.name(gammalib::strip_whitespace(par->attribute("name"))); + scale.has_grad(true); m_scales.push_back(scale); } diff --git a/src/model/GModelSky.cpp b/src/model/GModelSky.cpp index 8b8776a..3a067a8 100644 --- a/src/model/GModelSky.cpp +++ b/src/model/GModelSky.cpp @@ -1,7 +1,7 @@ /*************************************************************************** * GModelSky.cpp - Sky model class * * ----------------------------------------------------------------------- * - * copyright (C) 2011-2018 by Juergen Knoedlseder * + * copyright (C) 2011-2019 by Juergen Knoedlseder * * ----------------------------------------------------------------------- * * * * This program is free software: you can redistribute it and/or modify * @@ -1074,6 +1074,11 @@ std::string GModelSky::print(const GChatter& chatter) const for (int i = 0; i < n_temporal; ++i) { result.append("\n"+(*temporal())[i].print()); } + result.append("\n"+gammalib::parformat("Number of scale par's")); + result.append(gammalib::str(scales())); + for (int i = 0; i < scales(); ++i) { + result.append("\n"+scale(i).print()); + } } // endif: chatter was not silent @@ -1163,6 +1168,7 @@ void GModelSky::set_pointers(void) int n_spectral = (spectral() != NULL) ? spectral()->size() : 0; int n_temporal = (temporal() != NULL) ? temporal()->size() : 0; int n_pars = n_spatial + n_spectral + n_temporal; + int n_scales = m_scales.size(); // Continue only if there are parameters if (n_pars > 0) { @@ -1184,6 +1190,13 @@ void GModelSky::set_pointers(void) } + // Append scales if they exist + if (n_scales > 0) { + for (int i = 0; i < n_scales; ++i) { + m_pars.push_back(&m_scales[i]); + } + } + // Return return; } diff --git a/src/obs/GObservation.cpp b/src/obs/GObservation.cpp index 3ba9629..1035c4a 100644 --- a/src/obs/GObservation.cpp +++ b/src/obs/GObservation.cpp @@ -397,7 +397,7 @@ double GObservation::npred(const GModels& models, GVector* gradient) const #if defined(G_RANGE_CHECK) if (gradient != NULL) { if (models.npars() != gradient->size()) { - throw GException::gradient_par_mismatch(G_MODEL, + throw GException::gradient_par_mismatch(G_NPRED, gradient->size(), models.npars()); } diff --git a/src/obs/GResponse.cpp b/src/obs/GResponse.cpp index a669ceb..1f51031 100644 --- a/src/obs/GResponse.cpp +++ b/src/obs/GResponse.cpp @@ -1,7 +1,7 @@ /*************************************************************************** * GResponse.cpp - Abstract response base class * * ----------------------------------------------------------------------- * - * copyright (C) 2008-2018 by Juergen Knoedlseder * + * copyright (C) 2008-2019 by Juergen Knoedlseder * * ----------------------------------------------------------------------- * * * * This program is free software: you can redistribute it and/or modify * @@ -327,10 +327,9 @@ double GResponse::eval_prob(const GModelSky& model, // Continue only if IRF value is positive if (irf > 0.0) { - // If required, apply instrument specific model scaling - if (model.has_scales()) { - irf *= model.scale(obs.instrument()).value(); - } + // Optionally get scaling + double scale = (model.has_scales()) + ? model.scale(obs.instrument()).value() : 1.0; // Evaluate spectral and temporal components double spec = (model.spectral() != NULL) @@ -339,14 +338,14 @@ double GResponse::eval_prob(const GModelSky& model, ? model.temporal()->eval(srcTime, grad) : 1.0; // Compute probability - prob = spec * temp * irf; + prob = spec * temp * irf * scale; // Optionally compute partial derivatives if (grad) { // Multiply factors to spectral gradients if (model.spectral() != NULL) { - double fact = temp * irf; + double fact = temp * irf * scale; if (fact != 1.0) { for (int i = 0; i < model.spectral()->size(); ++i) { (*model.spectral())[i].factor_gradient((*model.spectral())[i].factor_gradient() * fact); @@ -356,7 +355,7 @@ double GResponse::eval_prob(const GModelSky& model, // Multiply factors to temporal gradients if (model.temporal() != NULL) { - double fact = spec * irf; + double fact = spec * irf * scale; if (fact != 1.0) { for (int i = 0; i < model.temporal()->size(); ++i) { (*model.temporal())[i].factor_gradient((*model.temporal())[i].factor_gradient() * fact); @@ -364,6 +363,19 @@ double GResponse::eval_prob(const GModelSky& model, } } + // Optionally compute scale gradient for instrument + if (model.has_scales()) { + for (int i = 0; i < model.scales(); ++i) { + double g_scale = 0.0; + if (model.scale(i).name() == obs.instrument()) { + if (model.scale(i).is_free()) { + g_scale = model.scale(i).scale() * spec * temp * irf; + } + } + model.scale(i).factor_gradient(g_scale); + } + } + } // endif: computed partial derivatives // Compile option: Check for NaN/Inf @@ -399,6 +411,13 @@ double GResponse::eval_prob(const GModelSky& model, } } + // Optionally reset scale gradients + if (model.has_scales()) { + for (int i = 0; i < model.scales(); ++i) { + model.scale(i).factor_gradient(0.0); + } + } + } // endelse: IRF value was not positive } // endif: Gamma-ray source model had a spatial component diff --git a/test/test_GModel.cpp b/test/test_GModel.cpp index bf7f21a..02eeecb 100644 --- a/test/test_GModel.cpp +++ b/test/test_GModel.cpp @@ -1,7 +1,7 @@ /*************************************************************************** * test_GModel.cpp - test GModel class * * ----------------------------------------------------------------------- * - * copyright (C) 2009-2017 by Juergen Knoedlseder * + * copyright (C) 2009-2019 by Juergen Knoedlseder * * ----------------------------------------------------------------------- * * * * This program is free software: you can redistribute it and/or modify * @@ -2474,7 +2474,7 @@ void TestGModel::test_xml_model(const std::string& name, void TestGModel::test_sky_model_content(const GModelSky& model) { // Test content of sky model - test_value(model.size(), 6, "Check number of model parameters"); + test_value(model.size(), 8, "Check number of model parameters"); test_value(model.name(), "1FGL J0005.7+3815", "Check source name"); test_value(model.instruments(), "LAT,CTA", "Check instruments"); test_value(model.ids(), "000001", "Check observation identifier"); @@ -2482,6 +2482,9 @@ void TestGModel::test_sky_model_content(const GModelSky& model) test_assert(model.tscalc(), "Check whether TS computation is requested"); test_assert(model.has_scales(), "Check whether sky model has scales"); test_assert(model.has_ts(), "Check whether sky model has TS value"); + test_value(model.scales(), 2, "Check number of instrument parameters"); + test_value(model.scale(0).value(), 1.1, 1.0e-6, "Check first scale factor"); + test_value(model.scale(1).value(), 0.5, 1.0e-6, "Check second scale factor"); test_value(model.scale("LAT").value(), 1.1, 1.0e-6, "Check LAT scale factor"); test_value(model.scale("CTA").value(), 0.5, 1.0e-6, "Check CTA scale factor"); test_value(model.type(), "PointSource", "Check source type"); -- libgit2 0.21.4