Nathan Kelley-Hoskins / gammalib

Commit d9f99bde0291386964b345817bc5eb16ff600fed

Authored by Jürgen Knödlseder
1 parent 59dffdf0

Allow fitting of instrument scale factors (#2892)

Inline Side-by-side
Showing 11 changed files with 167 additions and 47 deletions
ChangeLog
  View file @d9f99bd
1   -2019-04-26
  1 +2019-05-22
2 2  
3 3 * Version 1.6.0 released
4 4 ========================
5 5  
  6 + Allow fitting of instrument scale factors (#2892)
6 7 Add GCTAModelSpatialLookup class to handle CTA background lookups (#2878)
7 8 Implement energy-dependent Gaussian CTA spatial background model (#2868)
8 9 Take instrument scale factors in On/Off analysis into account (#2867)
... ...
1 1 New Features and Important Changes in GammaLib 1.6.0
2 2  
3   -26 April 2019
  3 +22 May 2019
4 4  
5 5  
6 6 1. Introduction
... ... @@ -73,6 +73,8 @@ The following methods have been added:
73 73 - GFitsTable::remove(int&)
74 74 - GFitsTable::remove(std::string&)
75 75 - GObservation::has_events()
  76 +- GModel::scales()
  77 +- GModel::scale(int&)
76 78 - GModelSpatialPointSource::type(std::string&)
77 79 - GModelSpectralPlaw::type(std::string&)
78 80 - GModelSpatial::has_free_pars()
... ... @@ -217,6 +219,9 @@ None
217 219  
218 220 8. Model module
219 221 ---------------
  222 +Add GModel::scales() method to return number of scales for this model and
  223 +add GModel::scale() parameter access methods (#2892).
  224 +
220 225 Add GModelSpatial::has_free_pars() method that checks whether a spatial
221 226 model component has free model parameters or not (#2715).
222 227  
... ...
doc/source/admin/release_history/1.6.rst
  View file @d9f99bd
... ... @@ -55,6 +55,8 @@ Bug fixes
55 55 Improvements
56 56 ------------
57 57  
  58 +* [`2892 <https://cta-redmine.irap.omp.eu/issues/2892>`_] -
  59 + Allow fitting of instrument scale factors
58 60 * [`2878 <https://cta-redmine.irap.omp.eu/issues/2878>`_] -
59 61 Add ``GCTAModelSpatialLookup`` class to handle CTA background lookups
60 62 * [`2868 <https://cta-redmine.irap.omp.eu/issues/2868>`_] -
... ...
include/GModel.hpp
  View file @d9f99bd
1 1 /***************************************************************************
2 2 * GModel.hpp - Abstract virtual model base class *
3 3 * ----------------------------------------------------------------------- *
4   - * copyright (C) 2009-2016 by Juergen Knoedlseder *
  4 + * copyright (C) 2009-2019 by Juergen Knoedlseder *
5 5 * ----------------------------------------------------------------------- *
6 6 * *
7 7 * This program is free software: you can redistribute it and/or modify *
... ... @@ -127,6 +127,7 @@ public:
127 127  
128 128 // Implemented methods
129 129 int size(void) const;
  130 + int scales(void) const;
130 131 GModelPar& at(const int& index);
131 132 const GModelPar& at(const int& index) const;
132 133 bool has_par(const std::string& name) const;
... ... @@ -140,6 +141,8 @@ public:
140 141 const bool& has_ts(void) const;
141 142 std::string instruments(void) const;
142 143 void instruments(const std::string& instruments);
  144 + GModelPar& scale(const int& index);
  145 + const GModelPar& scale(const int& index) const;
143 146 GModelPar scale(const std::string& instrument) const;
144 147 void scale(const GModelPar& par);
145 148 std::string ids(void) const;
... ... @@ -218,6 +221,20 @@ int GModel::size(void) const
218 221  
219 222  
220 223 /***********************************************************************//**
  224 + * @brief Return number of scale parameters in model
  225 + *
  226 + * @return Number of scale parameters in model.
  227 + *
  228 + * Returns the number of scale parameters in the model.
  229 + ***************************************************************************/
  230 +inline
  231 +int GModel::scales(void) const
  232 +{
  233 + return (int)m_scales.size();
  234 +}
  235 +
  236 +
  237 +/***********************************************************************//**
221 238 * @brief Return parameter name
222 239 *
223 240 * @return Parameter name.
... ...
inst/cta/src/GCTAOnOffObservation.cpp
  View file @d9f99bd
... ... @@ -2329,6 +2329,16 @@ double GCTAOnOffObservation::N_gamma(const GModels&amp; models,
2329 2329 double rmf_sum = 0.0;
2330 2330 #endif
2331 2331  
  2332 + // Get indices of spectral gradients. We need this to update
  2333 + // only the spectral parameter gradients later.
  2334 + std::vector<int> inx_spec;
  2335 + for (int k = 0; k < mptr->size(); ++k) {
  2336 + const GModelPar& par = (*mptr)[k];
  2337 + if (spectral->has_par(par.name())) {
  2338 + inx_spec.push_back(k);
  2339 + }
  2340 + }
  2341 +
2332 2342 // Set instrument scale factor
2333 2343 double scale = (sky->has_scales()) ? sky->scale(instrument()).value() : 1.0;
2334 2344  
... ... @@ -2358,16 +2368,20 @@ double GCTAOnOffObservation::N_gamma(const GModels&amp; models,
2358 2368 GEnergy etruemean = m_arf.ebounds().elogmean(itrue);
2359 2369 double etruewidth = m_arf.ebounds().ewidth(itrue).MeV();
2360 2370  
2361   - // Compute normalisation factors
2362   - double exposure = m_on_spec.exposure();
2363   - double norm_flux = arf * exposure * rmf * scale;
2364   - double norm_grad = norm_flux * etruewidth;
  2371 + // Compute normalisation factors. The gradient normalisation
  2372 + // factor uses the energy bin width, since the eval() method
  2373 + // returns a differential flux
  2374 + double exposure = m_on_spec.exposure();
  2375 + double norm_scale = arf * exposure * rmf;
  2376 + double norm_flux = norm_scale * scale;
  2377 + double norm_grad = norm_flux * etruewidth;
2365 2378  
2366 2379 // Determine number of gamma-ray events in model by
2367 2380 // computing the flux over the true energy bin in
2368 2381 // ph/cm2/s and multiplying this by effective area (cm2),
2369   - // livetime (s) and redistribution probability
2370   - value += spectral->flux(etruemin, etruemax) * norm_flux;
  2382 + // livetime (s), redistribution probability and scale
  2383 + double flux = spectral->flux(etruemin, etruemax);
  2384 + value += flux * norm_flux;
2371 2385  
2372 2386 // Determine the model gradients at the current true
2373 2387 // energy. The eval() method needs a time in case that the
... ... @@ -2375,11 +2389,22 @@ double GCTAOnOffObservation::N_gamma(const GModels&amp; models,
2375 2389 // dummy time here.
2376 2390 spectral->eval(etruemean, GTime(), true);
2377 2391  
2378   - // Compute the parameter gradients for all model parameters
2379   - for (int k = 0; k < mptr->size(); ++k) {
2380   - const GModelPar& par = (*mptr)[k];
2381   - if (par.is_free() && ipar+k < npars) {
2382   - (*grad)[ipar+k] += par.factor_gradient() * norm_grad;
  2392 + // Compute the parameter gradients for spectral parameters
  2393 + for (int k = 0; k < inx_spec.size(); ++k) {
  2394 + const GModelPar& par = (*mptr)[inx_spec[k]];
  2395 + if (par.is_free() && ipar+k < npars) {
  2396 + (*grad)[ipar+inx_spec[k]] += par.factor_gradient() * norm_grad;
  2397 + }
  2398 + }
  2399 +
  2400 + // Optionally compute scaling parameter gradient
  2401 + if (sky->has_scales()) {
  2402 + for (int k = 0; k < mptr->size(); ++k) {
  2403 + const GModelPar& par = (*mptr)[k];
  2404 + if ((par.name() == instrument()) && par.is_free() &&
  2405 + (ipar+k < npars)) {
  2406 + (*grad)[ipar+k] += par.scale() * flux * norm_scale;
  2407 + }
2383 2408 }
2384 2409 }
2385 2410  
... ...
pyext/GModel.i
  View file @d9f99bd
1 1 /***************************************************************************
2 2 * GModel.i - Abstract virtual model base class *
3 3 * ----------------------------------------------------------------------- *
4   - * copyright (C) 2009-2018 by Juergen Knoedlseder *
  4 + * copyright (C) 2009-2019 by Juergen Knoedlseder *
5 5 * ----------------------------------------------------------------------- *
6 6 * *
7 7 * This program is free software: you can redistribute it and/or modify *
... ... @@ -60,17 +60,19 @@ public:
60 60  
61 61 // Implemented methods
62 62 int size(void) const;
  63 + int scales(void) const;
63 64 bool has_par(const std::string& name) const;
64 65 bool has_scales(void) const;
65 66 const std::string& name(void) const;
66 67 void name(const std::string& name);
67   - std::string instruments(void) const;
68   - void instruments(const std::string& instruments);
69 68 const double& ts(void) const;
70 69 void ts(const double& ts);
71 70 const bool& tscalc(void) const;
72 71 void tscalc(const bool& tscalc);
73 72 const bool& has_ts(void) const;
  73 + std::string instruments(void) const;
  74 + void instruments(const std::string& instruments);
  75 + GModelPar& scale(const int& index);
74 76 GModelPar scale(const std::string& instrument) const;
75 77 void scale(const GModelPar& par);
76 78 std::string ids(void) const;
... ...
src/model/GModel.cpp
  View file @d9f99bd
1 1 /***************************************************************************
2 2 * GModel.cpp - Abstract virtual model base class *
3 3 * ----------------------------------------------------------------------- *
4   - * copyright (C) 2009-2015 by Juergen Knoedlseder *
  4 + * copyright (C) 2009-2019 by Juergen Knoedlseder *
5 5 * ----------------------------------------------------------------------- *
6 6 * *
7 7 * This program is free software: you can redistribute it and/or modify *
... ... @@ -35,6 +35,7 @@
35 35 /* __ Method name definitions ____________________________________________ */
36 36 #define G_ACCESS "GModel::operator[](std::string&)"
37 37 #define G_AT "GModel::at(int&)"
  38 +#define G_SCALE "GModelPar& GModel::scale(int&)"
38 39 #define G_WRITE_SCALES "GModel::write_scales(GXmlElement&)"
39 40  
40 41 /* __ Macros _____________________________________________________________ */
... ... @@ -348,6 +349,54 @@ void GModel::instruments(const std::string&amp; instruments)
348 349  
349 350  
350 351 /***********************************************************************//**
  352 + * @brief Returns reference to scale parameter by index
  353 + *
  354 + * @param[in] index Parameter index [0,...,size()-1].
  355 + * @return Reference to scale parameter.
  356 + *
  357 + * @exception GException::out_of_range
  358 + * Scale parameter index is out of range.
  359 + *
  360 + * Returns a reference to the scale parameter of the specified @p index.
  361 + * Throws an exception if @p index is not valid.
  362 + ***************************************************************************/
  363 +GModelPar& GModel::scale(const int& index)
  364 +{
  365 + // Raise exception if index is out of range
  366 + if (index < 0 || index >= scales()) {
  367 + throw GException::out_of_range(G_SCALE, index, scales());
  368 + }
  369 +
  370 + // Return reference
  371 + return (m_scales[index]);
  372 +}
  373 +
  374 +
  375 +/***********************************************************************//**
  376 + * @brief Returns reference to scale parameter by index (const version)
  377 + *
  378 + * @param[in] index Parameter index [0,...,size()-1].
  379 + * @return Reference to scale parameter.
  380 + *
  381 + * @exception GException::out_of_range
  382 + * Scale parameter index is out of range.
  383 + *
  384 + * Returns a reference to the scale parameter of the specified @p index.
  385 + * Throws an exception if @p index is not valid.
  386 + ***************************************************************************/
  387 +const GModelPar& GModel::scale(const int& index) const
  388 +{
  389 + // Raise exception if index is out of range
  390 + if (index < 0 || index >= scales()) {
  391 + throw GException::out_of_range(G_SCALE, index, scales());
  392 + }
  393 +
  394 + // Return reference
  395 + return (m_scales[index]);
  396 +}
  397 +
  398 +
  399 +/***********************************************************************//**
351 400 * @brief Returns model scale factor for a given instrument
352 401 *
353 402 * @param[in] instrument Instrument.
... ... @@ -715,23 +764,6 @@ std::string GModel::print_attributes(void) const
715 764 result.append("Computation requested");
716 765 }
717 766  
718   - // Append instrument scale factors
719   - result.append("\n"+gammalib::parformat("Instrument scale factors"));
720   - if (!m_scales.empty()) {
721   - for (int i = 0; i < m_scales.size(); ++i) {
722   - if (i > 0) {
723   - result.append(", ");
724   - }
725   - result.append(m_scales[i].name());
726   - result.append("=");
727   - result.append(gammalib::str(m_scales[i].value()));
728   - }
729   - result.append(", others unity");
730   - }
731   - else {
732   - result.append("unity");
733   - }
734   -
735 767 // Append observation identifiers
736 768 result.append("\n"+gammalib::parformat("Observation identifiers"));
737 769 if (!m_ids.empty()) {
... ... @@ -788,6 +820,7 @@ void GModel::read_scales(const GXmlElement&amp; xml)
788 820 GModelPar scale;
789 821 scale.read(*par);
790 822 scale.name(gammalib::strip_whitespace(par->attribute("name")));
  823 + scale.has_grad(true);
791 824 m_scales.push_back(scale);
792 825 }
793 826  
... ...
src/model/GModelSky.cpp
  View file @d9f99bd
1 1 /***************************************************************************
2 2 * GModelSky.cpp - Sky model class *
3 3 * ----------------------------------------------------------------------- *
4   - * copyright (C) 2011-2018 by Juergen Knoedlseder *
  4 + * copyright (C) 2011-2019 by Juergen Knoedlseder *
5 5 * ----------------------------------------------------------------------- *
6 6 * *
7 7 * This program is free software: you can redistribute it and/or modify *
... ... @@ -1074,6 +1074,11 @@ std::string GModelSky::print(const GChatter&amp; chatter) const
1074 1074 for (int i = 0; i < n_temporal; ++i) {
1075 1075 result.append("\n"+(*temporal())[i].print());
1076 1076 }
  1077 + result.append("\n"+gammalib::parformat("Number of scale par's"));
  1078 + result.append(gammalib::str(scales()));
  1079 + for (int i = 0; i < scales(); ++i) {
  1080 + result.append("\n"+scale(i).print());
  1081 + }
1077 1082  
1078 1083 } // endif: chatter was not silent
1079 1084  
... ... @@ -1163,6 +1168,7 @@ void GModelSky::set_pointers(void)
1163 1168 int n_spectral = (spectral() != NULL) ? spectral()->size() : 0;
1164 1169 int n_temporal = (temporal() != NULL) ? temporal()->size() : 0;
1165 1170 int n_pars = n_spatial + n_spectral + n_temporal;
  1171 + int n_scales = m_scales.size();
1166 1172  
1167 1173 // Continue only if there are parameters
1168 1174 if (n_pars > 0) {
... ... @@ -1184,6 +1190,13 @@ void GModelSky::set_pointers(void)
1184 1190  
1185 1191 }
1186 1192  
  1193 + // Append scales if they exist
  1194 + if (n_scales > 0) {
  1195 + for (int i = 0; i < n_scales; ++i) {
  1196 + m_pars.push_back(&m_scales[i]);
  1197 + }
  1198 + }
  1199 +
1187 1200 // Return
1188 1201 return;
1189 1202 }
... ...
src/obs/GObservation.cpp
  View file @d9f99bd
... ... @@ -397,7 +397,7 @@ double GObservation::npred(const GModels&amp; models, GVector* gradient) const
397 397 #if defined(G_RANGE_CHECK)
398 398 if (gradient != NULL) {
399 399 if (models.npars() != gradient->size()) {
400   - throw GException::gradient_par_mismatch(G_MODEL,
  400 + throw GException::gradient_par_mismatch(G_NPRED,
401 401 gradient->size(),
402 402 models.npars());
403 403 }
... ...
src/obs/GResponse.cpp
  View file @d9f99bd
1 1 /***************************************************************************
2 2 * GResponse.cpp - Abstract response base class *
3 3 * ----------------------------------------------------------------------- *
4   - * copyright (C) 2008-2018 by Juergen Knoedlseder *
  4 + * copyright (C) 2008-2019 by Juergen Knoedlseder *
5 5 * ----------------------------------------------------------------------- *
6 6 * *
7 7 * This program is free software: you can redistribute it and/or modify *
... ... @@ -327,10 +327,9 @@ double GResponse::eval_prob(const GModelSky&amp; model,
327 327 // Continue only if IRF value is positive
328 328 if (irf > 0.0) {
329 329  
330   - // If required, apply instrument specific model scaling
331   - if (model.has_scales()) {
332   - irf *= model.scale(obs.instrument()).value();
333   - }
  330 + // Optionally get scaling
  331 + double scale = (model.has_scales())
  332 + ? model.scale(obs.instrument()).value() : 1.0;
334 333  
335 334 // Evaluate spectral and temporal components
336 335 double spec = (model.spectral() != NULL)
... ... @@ -339,14 +338,14 @@ double GResponse::eval_prob(const GModelSky&amp; model,
339 338 ? model.temporal()->eval(srcTime, grad) : 1.0;
340 339  
341 340 // Compute probability
342   - prob = spec * temp * irf;
  341 + prob = spec * temp * irf * scale;
343 342  
344 343 // Optionally compute partial derivatives
345 344 if (grad) {
346 345  
347 346 // Multiply factors to spectral gradients
348 347 if (model.spectral() != NULL) {
349   - double fact = temp * irf;
  348 + double fact = temp * irf * scale;
350 349 if (fact != 1.0) {
351 350 for (int i = 0; i < model.spectral()->size(); ++i) {
352 351 (*model.spectral())[i].factor_gradient((*model.spectral())[i].factor_gradient() * fact);
... ... @@ -356,7 +355,7 @@ double GResponse::eval_prob(const GModelSky&amp; model,
356 355  
357 356 // Multiply factors to temporal gradients
358 357 if (model.temporal() != NULL) {
359   - double fact = spec * irf;
  358 + double fact = spec * irf * scale;
360 359 if (fact != 1.0) {
361 360 for (int i = 0; i < model.temporal()->size(); ++i) {
362 361 (*model.temporal())[i].factor_gradient((*model.temporal())[i].factor_gradient() * fact);
... ... @@ -364,6 +363,19 @@ double GResponse::eval_prob(const GModelSky&amp; model,
364 363 }
365 364 }
366 365  
  366 + // Optionally compute scale gradient for instrument
  367 + if (model.has_scales()) {
  368 + for (int i = 0; i < model.scales(); ++i) {
  369 + double g_scale = 0.0;
  370 + if (model.scale(i).name() == obs.instrument()) {
  371 + if (model.scale(i).is_free()) {
  372 + g_scale = model.scale(i).scale() * spec * temp * irf;
  373 + }
  374 + }
  375 + model.scale(i).factor_gradient(g_scale);
  376 + }
  377 + }
  378 +
367 379 } // endif: computed partial derivatives
368 380  
369 381 // Compile option: Check for NaN/Inf
... ... @@ -399,6 +411,13 @@ double GResponse::eval_prob(const GModelSky&amp; model,
399 411 }
400 412 }
401 413  
  414 + // Optionally reset scale gradients
  415 + if (model.has_scales()) {
  416 + for (int i = 0; i < model.scales(); ++i) {
  417 + model.scale(i).factor_gradient(0.0);
  418 + }
  419 + }
  420 +
402 421 } // endelse: IRF value was not positive
403 422  
404 423 } // endif: Gamma-ray source model had a spatial component
... ...
test/test_GModel.cpp
  View file @d9f99bd
1 1 /***************************************************************************
2 2 * test_GModel.cpp - test GModel class *
3 3 * ----------------------------------------------------------------------- *
4   - * copyright (C) 2009-2017 by Juergen Knoedlseder *
  4 + * copyright (C) 2009-2019 by Juergen Knoedlseder *
5 5 * ----------------------------------------------------------------------- *
6 6 * *
7 7 * This program is free software: you can redistribute it and/or modify *
... ... @@ -2474,7 +2474,7 @@ void TestGModel::test_xml_model(const std::string&amp; name,
2474 2474 void TestGModel::test_sky_model_content(const GModelSky& model)
2475 2475 {
2476 2476 // Test content of sky model
2477   - test_value(model.size(), 6, "Check number of model parameters");
  2477 + test_value(model.size(), 8, "Check number of model parameters");
2478 2478 test_value(model.name(), "1FGL J0005.7+3815", "Check source name");
2479 2479 test_value(model.instruments(), "LAT,CTA", "Check instruments");
2480 2480 test_value(model.ids(), "000001", "Check observation identifier");
... ... @@ -2482,6 +2482,9 @@ void TestGModel::test_sky_model_content(const GModelSky&amp; model)
2482 2482 test_assert(model.tscalc(), "Check whether TS computation is requested");
2483 2483 test_assert(model.has_scales(), "Check whether sky model has scales");
2484 2484 test_assert(model.has_ts(), "Check whether sky model has TS value");
  2485 + test_value(model.scales(), 2, "Check number of instrument parameters");
  2486 + test_value(model.scale(0).value(), 1.1, 1.0e-6, "Check first scale factor");
  2487 + test_value(model.scale(1).value(), 0.5, 1.0e-6, "Check second scale factor");
2485 2488 test_value(model.scale("LAT").value(), 1.1, 1.0e-6, "Check LAT scale factor");
2486 2489 test_value(model.scale("CTA").value(), 0.5, 1.0e-6, "Check CTA scale factor");
2487 2490 test_value(model.type(), "PointSource", "Check source type");
... ...