Francois Brun / gammalib

Commit eb6dc1336ca5213c509b94329f54b57ff3509b7d

Authored by Jürgen Knödlseder
1 parent 2c5c39b5

Integrate multiplicative spectral model

Inline Side-by-side
Showing 6 changed files with 366 additions and 311 deletions
ChangeLog
  View file @eb6dc13
1   -2016-10-05
  1 +2016-10-06
2 2  
3 3 * Version 1.2.0 released
4 4 =========================
5 5  
  6 + Add multiplicative spectral model (#1861)
6 7 Add additional time setting and getting methods (#1856)
7 8 Check CTA response table dimension for consistency (#1848)
8 9 Remove eval_gradients() methods and remplace by eval() method with optional gradients argument (#1847)
... ...
1 1 New Features and Important Changes in GammaLib 1.2.0
2 2  
3   -5 October 2016
  3 +6 October 2016
4 4  
5 5  
6 6 1. Introduction
... ... @@ -14,6 +14,7 @@ made since the last release of GammaLib.
14 14 The following classes have been added:
15 15 - GModelSpectralExpInvPlaw
16 16 - GModelSpectralPlawEnergyFlux
  17 +- GModelSpectralMultiplicative
17 18  
18 19 The following classes have been removed:
19 20 - ...
... ... @@ -143,6 +144,9 @@ parameter instead of a cut-off energy for the exponentially cut-off
143 144 power law. The cut-off parameter is the inverse of the cut-off energy
144 145 (#1800).
145 146  
  147 +The GModelSpectralMultiplicative class has been added that implements the
  148 +product of spectral model components as a new spectral model (#1861).
  149 +
146 150  
147 151 9. Numerics module
148 152 ------------------
... ...
include/GModelSpectralMultiplicative.hpp
  View file @eb6dc13
1 1 /***************************************************************************
2   - * GModelSpectralMultiplicative.hpp - Spectral power law model class *
  2 + * GModelSpectralMultiplicative.hpp - Multiplicative spectral model class *
3 3 * ----------------------------------------------------------------------- *
4   - * copyright (C) 2016 by Michael Mayer *
  4 + * copyright (C) 2016 by Michael Mayer *
5 5 * ----------------------------------------------------------------------- *
6 6 * *
7 7 * This program is free software: you can redistribute it and/or modify *
... ... @@ -41,16 +41,14 @@ class GTime;
41 41 class GXmlElement;
42 42  
43 43  
44   -
45 44 /***********************************************************************//**
46 45 * @class GModelSpectralMultiplicative
47 46 *
48 47 * @brief Multiplicative spectral model class
49 48 *
50   - * This class implements a Multiplicative spectrum. The model is defined by the sum
51   - * of several individual spectral models. Each spectral model can be added via
52   - * the XML interface or using the append() method.
53   - *
  49 + * This class implements a spectral model that is a multiplication spectral
  50 + * model components. The spectral model components can be defined in an
  51 + * XML file, or added using the append() method.
54 52 ***************************************************************************/
55 53 class GModelSpectralMultiplicative : public GModelSpectral {
56 54  
... ... @@ -65,24 +63,24 @@ public:
65 63 virtual GModelSpectralMultiplicative& operator=(const GModelSpectralMultiplicative& model);
66 64  
67 65 // Implemented pure virtual methods
68   - virtual void clear(void);
  66 + virtual void clear(void);
69 67 virtual GModelSpectralMultiplicative* clone(void) const;
70   - virtual std::string classname(void) const;
71   - virtual std::string type(void) const;
72   - virtual double eval(const GEnergy& srcEng,
73   - const GTime& srcTime = GTime(),
74   - const bool& gradients = false) const;
75   - virtual double flux(const GEnergy& emin,
76   - const GEnergy& emax) const;
77   - virtual double eflux(const GEnergy& emin,
78   - const GEnergy& emax) const;
79   - virtual GEnergy mc(const GEnergy& emin,
80   - const GEnergy& emax,
81   - const GTime& time,
82   - GRan& ran) const;
83   - virtual void read(const GXmlElement& xml);
84   - virtual void write(GXmlElement& xml) const;
85   - virtual std::string print(const GChatter& chatter = NORMAL) const;
  68 + virtual std::string classname(void) const;
  69 + virtual std::string type(void) const;
  70 + virtual double eval(const GEnergy& srcEng,
  71 + const GTime& srcTime = GTime(),
  72 + const bool& gradients = false) const;
  73 + virtual double flux(const GEnergy& emin,
  74 + const GEnergy& emax) const;
  75 + virtual double eflux(const GEnergy& emin,
  76 + const GEnergy& emax) const;
  77 + virtual GEnergy mc(const GEnergy& emin,
  78 + const GEnergy& emax,
  79 + const GTime& time,
  80 + GRan& ran) const;
  81 + virtual void read(const GXmlElement& xml);
  82 + virtual void write(GXmlElement& xml) const;
  83 + virtual std::string print(const GChatter& chatter = NORMAL) const;
86 84  
87 85 // Other methods
88 86 void append(const GModelSpectral& spec, const std::string& name="");
... ... @@ -98,34 +96,32 @@ protected:
98 96 void add_component(const GXmlElement& spec);
99 97 void update_mc_cache(const GEnergy& emin, const GEnergy& emax) const;
100 98  
101   - // Class to determine to the integral photon flux
102   - class flux_kern : public GFunction {
103   - public:
104   - // Constructor
105   - flux_kern(std::vector<GModelSpectral*> spec) :
106   - m_spec(spec){};
107   -
108   - // Method
109   - double eval(const double& x) {
110   - GEnergy energy(x, "MeV");
111   - double value = m_spec[0]->eval(energy);
112   - for(int i=1;i<m_spec.size();++i){
113   - value*=m_spec[i]->eval(energy);
114   - }
115   - return value;
116   - }
117   -
118   - // Members
119   - protected:
120   - std::vector<GModelSpectral*> m_spec; //!< Spectral models
121   - };
122   -
123   - // Class to determine the integral energy flux, derviation of flux_kern
  99 + // Class to determine the integral photon flux
  100 + class flux_kern : public GFunction {
  101 + public:
  102 + // Constructor
  103 + flux_kern(std::vector<GModelSpectral*> spec) : m_spec(spec) {}
  104 +
  105 + // Method
  106 + double eval(const double& x) {
  107 + GEnergy energy(x, "MeV");
  108 + double value = m_spec[0]->eval(energy);
  109 + for (int i = 1 ; i < m_spec.size(); ++i) {
  110 + value *= m_spec[i]->eval(energy);
  111 + }
  112 + return value;
  113 + }
  114 +
  115 + protected:
  116 + // Protected members
  117 + std::vector<GModelSpectral*> m_spec; //!< Spectral models
  118 + };
  119 +
  120 + // Class to determine the integral energy flux, derived from flux_kern
124 121 class eflux_kern : public flux_kern {
125 122 public:
126 123 // Constructor
127   - eflux_kern(std::vector<GModelSpectral*> spec) :
128   - flux_kern(spec) {};
  124 + eflux_kern(std::vector<GModelSpectral*> spec) : flux_kern(spec) {}
129 125  
130 126 // Method
131 127 double eval(const double& x) {
... ... @@ -140,10 +136,9 @@ protected:
140 136  
141 137 // MC cache
142 138 mutable GModelSpectralNodes m_mc_spectrum; //!< MC spectrum cache
143   - mutable int m_n_mc_nodes; //!< Probailities of individual components
144 139 mutable GEnergy m_mc_emin; //!< Last minimum energy
145 140 mutable GEnergy m_mc_emax; //!< Last maximum energy
146   -
  141 + mutable std::vector<double> m_mc_values; //!< Parameter values
147 142 };
148 143  
149 144  
... ... @@ -186,5 +181,4 @@ int GModelSpectralMultiplicative::components(void) const
186 181 return (m_spectral.size());
187 182 }
188 183  
189   -
190 184 #endif /* GMODELSPECTRALMULTIPLICATIVE_HPP */
... ...
pyext/GModelSpectralMultiplicative.i
  View file @eb6dc13
1 1 /***************************************************************************
2   - * GModelSpectralMultiplicative.i - Spectral power law model class *
  2 + * GModelSpectralMultiplicative.i - Multiplicative spectral model class *
3 3 * ----------------------------------------------------------------------- *
4 4 * copyright (C) 2016 Michael Mayer *
5 5 * ----------------------------------------------------------------------- *
... ... @@ -19,7 +19,8 @@
19 19 * *
20 20 ***************************************************************************/
21 21 /**
22   - * @file GModelSpecComposite spectral model class interface definition
  22 + * @file GModelSpectralMultiplicative.i
  23 + * @brief Multiplicative spectral model class interface definition
23 24 * @author Michael Mayer
24 25 */
25 26 %{
... ... @@ -31,7 +32,7 @@
31 32 /***********************************************************************//**
32 33 * @class GModelSpectralMultiplicative
33 34 *
34   - * @brief Composite spectral model class
  35 + * @brief Multiplicative spectral model class
35 36 ***************************************************************************/
36 37 class GModelSpectralMultiplicative : public GModelSpectral {
37 38 public:
... ... @@ -42,30 +43,29 @@ public:
42 43 virtual ~GModelSpectralMultiplicative(void);
43 44  
44 45 // Implemented pure virtual methods
45   - virtual void clear(void);
  46 + virtual void clear(void);
46 47 virtual GModelSpectralMultiplicative* clone(void) const;
47   - virtual std::string classname(void) const;
48   - virtual std::string type(void) const;
49   - virtual double eval(const GEnergy& srcEng,
50   - const GTime& srcTime = GTime(),
51   - const bool& gradients = false) const;
52   - virtual double flux(const GEnergy& emin,
53   - const GEnergy& emax) const;
54   - virtual double eflux(const GEnergy& emin,
55   - const GEnergy& emax) const;
56   - virtual GEnergy mc(const GEnergy& emin,
57   - const GEnergy& emax,
58   - const GTime& time,
59   - GRan& ran) const;
60   - virtual void read(const GXmlElement& xml);
61   - virtual void write(GXmlElement& xml) const;
  48 + virtual std::string classname(void) const;
  49 + virtual std::string type(void) const;
  50 + virtual double eval(const GEnergy& srcEng,
  51 + const GTime& srcTime = GTime(),
  52 + const bool& gradients = false) const;
  53 + virtual double flux(const GEnergy& emin,
  54 + const GEnergy& emax) const;
  55 + virtual double eflux(const GEnergy& emin,
  56 + const GEnergy& emax) const;
  57 + virtual GEnergy mc(const GEnergy& emin,
  58 + const GEnergy& emax,
  59 + const GTime& time,
  60 + GRan& ran) const;
  61 + virtual void read(const GXmlElement& xml);
  62 + virtual void write(GXmlElement& xml) const;
62 63  
63 64 // Other methods
64   - void append(const GModelSpectral& spec, const std::string& name="");
65   - int components(void) const;
66   - GModelSpectral* component(const int& index) const;
67   - GModelSpectral* component(const std::string& name) const;
68   -
  65 + void append(const GModelSpectral& spec, const std::string& name="");
  66 + int components(void) const;
  67 + const GModelSpectral* component(const int& index) const;
  68 + const GModelSpectral* component(const std::string& name) const;
69 69 };
70 70  
71 71  
... ...
src/model/GModelSpectralMultiplicative.cpp
  View file @eb6dc13
1 1 /***************************************************************************
2   - * GModelSpectralMultiplicative.cpp - Spectral multiplicative model class *
  2 + * GModelSpectralMultiplicative.cpp - Multiplicative spectral model class *
3 3 * ----------------------------------------------------------------------- *
4 4 * copyright (C) 2016 by Michael Mayer *
5 5 * ----------------------------------------------------------------------- *
... ... @@ -41,14 +41,16 @@
41 41 /* __ Constants __________________________________________________________ */
42 42  
43 43 /* __ Globals ____________________________________________________________ */
44   -const GModelSpectralMultiplicative g_spectral_multi_seed;
45   -const GModelSpectralRegistry g_spectral_multi_registry(&g_spectral_multi_seed);
  44 +const GModelSpectralMultiplicative g_spectral_multi_seed;
  45 +const GModelSpectralRegistry g_spectral_multi_registry(&g_spectral_multi_seed);
46 46  
47 47 /* __ Method name definitions ____________________________________________ */
48   -#define G_MC "GModelSpectralMultiplicative::mc(GEnergy&, GEnergy&, GTime&, GRan&)"
49   -#define G_WRITE "GModelSpectralMultiplicative::write(GXmlElement&)"
50   -#define G_COMPONENT_INDEX "GModelSpectralMultiplicative::component(const int&)"
51   -#define G_COMPONENT_NAME "GModelSpectralMultiplicative::component(const std::string&)"
  48 +#define G_MC "GModelSpectralMultiplicative::mc(GEnergy&, GEnergy&, GTime&, "\
  49 + "GRan&)"
  50 +#define G_WRITE "GModelSpectralMultiplicative::write(GXmlElement&)"
  51 +#define G_COMPONENT_INDEX "GModelSpectralMultiplicative::component(int&)"
  52 +#define G_COMPONENT_NAME "GModelSpectralMultiplicative::component"\
  53 + "(std::string&)"
52 54  
53 55 /* __ Macros _____________________________________________________________ */
54 56  
... ... @@ -66,7 +68,8 @@ const GModelSpectralRegistry g_spectral_multi_registry(&amp;g_spectral_multi_seed);
66 68 /***********************************************************************//**
67 69 * @brief Void constructor
68 70 ***************************************************************************/
69   -GModelSpectralMultiplicative::GModelSpectralMultiplicative(void) : GModelSpectral()
  71 +GModelSpectralMultiplicative::GModelSpectralMultiplicative(void) :
  72 + GModelSpectral()
70 73 {
71 74 // Initialise private members for clean destruction
72 75 init_members();
... ... @@ -79,13 +82,14 @@ GModelSpectralMultiplicative::GModelSpectralMultiplicative(void) : GModelSpectra
79 82 /***********************************************************************//**
80 83 * @brief XML constructor
81 84 *
82   - * @param[in] xml XML element containing position information.
  85 + * @param[in] xml XML element containing spectral model information.
83 86 *
84   - * Constructs a multiplicative spectral model by extracting information from an
85   - * XML element. See the read() method for more information about the expected
86   - * structure of the XML element.
  87 + * Constructs a multiplicative spectral model by extracting information from
  88 + * an XML element. See the read() method for more information about the
  89 + * expected structure of the XML element.
87 90 ***************************************************************************/
88   -GModelSpectralMultiplicative::GModelSpectralMultiplicative(const GXmlElement& xml) : GModelSpectral()
  91 +GModelSpectralMultiplicative::GModelSpectralMultiplicative(const GXmlElement& xml) :
  92 + GModelSpectral()
89 93 {
90 94 // Initialise members
91 95 init_members();
... ... @@ -103,8 +107,8 @@ GModelSpectralMultiplicative::GModelSpectralMultiplicative(const GXmlElement&amp; xm
103 107 *
104 108 * @param[in] model Multiplicative spectral model.
105 109 ***************************************************************************/
106   -GModelSpectralMultiplicative::GModelSpectralMultiplicative(const GModelSpectralMultiplicative& model)
107   - : GModelSpectral(model)
  110 +GModelSpectralMultiplicative::GModelSpectralMultiplicative(const GModelSpectralMultiplicative& model) :
  111 + GModelSpectral(model)
108 112 {
109 113 // Initialise members
110 114 init_members();
... ... @@ -139,8 +143,8 @@ GModelSpectralMultiplicative::~GModelSpectralMultiplicative(void)
139 143 /***********************************************************************//**
140 144 * @brief Assignment operator
141 145 *
142   - * @param[in] model Spectral power law model.
143   - * @return Spectral power law model.
  146 + * @param[in] model Multiplicative spectral model.
  147 + * @return Multiplicative spectral model.
144 148 ***************************************************************************/
145 149 GModelSpectralMultiplicative& GModelSpectralMultiplicative::operator=(const GModelSpectralMultiplicative& model)
146 150 {
... ... @@ -216,61 +220,65 @@ GModelSpectralMultiplicative* GModelSpectralMultiplicative::clone(void) const
216 220 * \prod_{i=0}^{N} {M_{\rm i}}(\rm srcEng, srcTime)
217 221 * \f]
218 222 *
219   - * where \f${M_{\rm i}}\f$ is the model evaluation of the i-th model
  223 + * where \f${M_{\rm i}}\f$ is the i-th model component.
  224 + *
  225 + * If the @p gradients flag is true the method will also compute the partial
  226 + * derivatives of each model component with respect to the parameters using
220 227 *
221   - * If the @p gradients flag is true the method will also compute the
222   - * partial derivatives of each model component with respect to the parameters using
  228 + * @todo Add formula for gradient computation
223 229 *
224 230 * @todo The method expects that energy!=0. Otherwise Inf or NaN may result.
225 231 ***************************************************************************/
226 232 double GModelSpectralMultiplicative::eval(const GEnergy& srcEng,
227   - const GTime& srcTime,
228   - const bool& gradients) const
  233 + const GTime& srcTime,
  234 + const bool& gradients) const
229 235 {
230 236 // Initialise result
231   - double value = 0.0;
  237 + double value = 0.0;
232 238  
233   - // Check for available spectral components
234   - if (m_spectral.size()) {
  239 + // Check for available spectral components
  240 + if (m_spectral.size()) {
235 241  
236   - // Add first spectral component
237   - value += m_spectral[0]->eval(srcEng, srcTime, gradients);
  242 + // Set first spectral component
  243 + value = m_spectral[0]->eval(srcEng, srcTime, gradients);
238 244  
239   - // Loop over model components
240   - for (int i = 1; i < m_spectral.size(); ++i) {
241   - value *= m_spectral[i]->eval(srcEng, srcTime, gradients);
242   - }
  245 + // Loop over model components
  246 + for (int i = 1; i < m_spectral.size(); ++i) {
  247 + value *= m_spectral[i]->eval(srcEng, srcTime, gradients);
  248 + }
  249 +
  250 + } // endfor: loop over model components
243 251  
244   - } // endfor: loop over model components
  252 + // Modify gradients if requested
  253 + if (gradients) {
245 254  
246   - // Modify gradients if requested
247   - if (gradients) {
  255 + // Loop over model components
  256 + for (int i = 0; i < m_spectral.size(); ++i) {
248 257  
249   - // Loop over models
250   - for(int i = 0; i < m_spectral.size(); ++i) {
  258 + // Initialise scaling factor
  259 + double factor = 1.0;
251 260  
252   - // Initialise scaling factor
253   - double factor = 1.0;
  261 + // Loop over other model components and compute factor
  262 + for (int j = 0; j < m_spectral.size(); ++j) {
  263 + if (i != j) {
  264 + factor *= m_spectral[j]->eval(srcEng, srcTime, false);
  265 + }
  266 + }
254 267  
255   - //Loop over other models and compute factor
256   - for(int j = 0; j < m_spectral.size(); ++j) {
257   - if (i != j) {
258   - factor *= m_spectral[j]->eval(srcEng, srcTime, false);
259   - }
260   - }
  268 + // Loop over model parameters
  269 + for (int ipar = 0; ipar < m_spectral[i]->size(); ++ipar) {
261 270  
262   - // Loop over model parameters
263   - for(int ipar = 0; ipar < m_spectral[i]->size(); ++ipar) {
  271 + // Get reference to model parameter
  272 + GModelPar& par = m_spectral[i]->operator[](ipar);
264 273  
265   - // Get reference to model parameter
266   - GModelPar& par = m_spectral[i]->operator[](ipar);
  274 + // Scale parameter gradient
  275 + par.gradient(par.gradient()*factor);
267 276  
268   - // Scale parameter gradient
269   - par.gradient(par.gradient()*factor);
  277 + } // endfor: loop over model parameters
270 278  
271   - } // endfor: loop over model parameters
272   - } // endfor: loop over models
273   - } //endif: compute gradients
  279 + } // endfor: loop over models
  280 +
  281 + } //endif: compute gradients
274 282  
275 283  
276 284 // Compile option: Check for NaN/Inf
... ... @@ -300,7 +308,7 @@ double GModelSpectralMultiplicative::eval(const GEnergy&amp; srcEng,
300 308 * Computes the photon flux of multiplicative spectral model
301 309 ***************************************************************************/
302 310 double GModelSpectralMultiplicative::flux(const GEnergy& emin,
303   - const GEnergy& emax) const
  311 + const GEnergy& emax) const
304 312 {
305 313 // Initialise flux
306 314 double flux = 0.0;
... ... @@ -308,17 +316,17 @@ double GModelSpectralMultiplicative::flux(const GEnergy&amp; emin,
308 316 // Compute only if integration range is valid and components are available
309 317 if (emin < emax && m_spectral.size()) {
310 318  
311   - // Initialise function to integrate
312   - flux_kern kernel(m_spectral);
  319 + // Initialise function to integrate
  320 + flux_kern kernel(m_spectral);
313 321  
314   - // Initialise integral class with function
315   - GIntegral integral(&kernel);
  322 + // Initialise integral class with function
  323 + GIntegral integral(&kernel);
316 324  
317   - // Set integration precision
318   - integral.eps(1.0e-8);
  325 + // Set integration precision
  326 + integral.eps(1.0e-8);
319 327  
320   - // Calculate integral between emin and emax
321   - flux = integral.romberg(emin.MeV(), emax.MeV());
  328 + // Calculate integral between emin and emax
  329 + flux = integral.romberg(emin.MeV(), emax.MeV());
322 330  
323 331 } // endif: integration range was valid
324 332  
... ... @@ -337,30 +345,30 @@ double GModelSpectralMultiplicative::flux(const GEnergy&amp; emin,
337 345 * Computes the energy flux of multiplicative spectral model
338 346 ***************************************************************************/
339 347 double GModelSpectralMultiplicative::eflux(const GEnergy& emin,
340   - const GEnergy& emax) const
  348 + const GEnergy& emax) const
341 349 {
342   - // Initialise eflux
343   - double eflux = 0.0;
  350 + // Initialise eflux
  351 + double eflux = 0.0;
344 352  
345   - // Compute only if integration range is valid and components are available
346   - if (emin < emax && m_spectral.size()) {
  353 + // Compute only if integration range is valid and components are available
  354 + if (emin < emax && m_spectral.size()) {
347 355  
348   - // Initialise function to integrate
349   - eflux_kern kernel(m_spectral);
  356 + // Initialise function to integrate
  357 + eflux_kern kernel(m_spectral);
350 358  
351   - // Initialise integral class with function
352   - GIntegral integral(&kernel);
  359 + // Initialise integral class with function
  360 + GIntegral integral(&kernel);
353 361  
354   - // Set integration precision
355   - integral.eps(1.0e-8);
  362 + // Set integration precision
  363 + integral.eps(1.0e-8);
356 364  
357   - // Calculate integral between emin and emax
358   - eflux = integral.romberg(emin.MeV(), emax.MeV());
  365 + // Calculate integral between emin and emax
  366 + eflux = integral.romberg(emin.MeV(), emax.MeV());
359 367  
360   - } // endif: integration range was valid
  368 + } // endif: integration range was valid
361 369  
362   - // Return flux
363   - return eflux;
  370 + // Return flux
  371 + return eflux;
364 372 }
365 373  
366 374  
... ... @@ -376,12 +384,13 @@ double GModelSpectralMultiplicative::eflux(const GEnergy&amp; emin,
376 384 * @exception GException::erange_invalid
377 385 * Energy range is invalid (emin < emax required).
378 386 *
379   - * Returns Monte Carlo energy by randomly drawing from a multiplicative spectral model.
  387 + * Returns Monte Carlo energy by randomly drawing from a multiplicative
  388 + * spectral model.
380 389 ***************************************************************************/
381 390 GEnergy GModelSpectralMultiplicative::mc(const GEnergy& emin,
382   - const GEnergy& emax,
383   - const GTime& time,
384   - GRan& ran) const
  391 + const GEnergy& emax,
  392 + const GTime& time,
  393 + GRan& ran) const
385 394 {
386 395 // Throw an exception if energy range is invalid
387 396 if (emin >= emax) {
... ... @@ -391,9 +400,9 @@ GEnergy GModelSpectralMultiplicative::mc(const GEnergy&amp; emin,
391 400  
392 401 // Throw exception if model container is empty
393 402 if (m_spectral.size() == 0) {
394   - throw GException::runtime_error(G_MC,
395   - "Multiplicative spectral model is empty. It is required to have at"\
396   - " least one spectral model to run simulations");
  403 + std::string msg = "Multiplicative spectral model is empty. At least"
  404 + " one spectral model is required for simulations.";
  405 + throw GException::runtime_error(G_MC, msg);
397 406 }
398 407  
399 408 // Update MC cache
... ... @@ -416,19 +425,19 @@ GEnergy GModelSpectralMultiplicative::mc(const GEnergy&amp; emin,
416 425 ***************************************************************************/
417 426 void GModelSpectralMultiplicative::read(const GXmlElement& xml)
418 427 {
419   - // Get number of spectral components
420   - int n_spectrals = xml.elements("spectrum");
  428 + // Get number of spectral components
  429 + int n_spectrals = xml.elements("spectrum");
421 430  
422   - // Loop over spectral elements
423   - for(int i = 0; i < n_spectrals; ++i) {
  431 + // Loop over spectral elements
  432 + for (int i = 0; i < n_spectrals; ++i) {
424 433  
425   - // Get spectral XML element
426   - const GXmlElement* spec = xml.element("spectrum", i);
  434 + // Get spectral XML element
  435 + const GXmlElement* spec = xml.element("spectrum", i);
427 436  
428   - // Add spectral component
429   - add_component(*spec);
  437 + // Add spectral component
  438 + add_component(*spec);
430 439  
431   - } // endfor: loop over components
  440 + } // endfor: loop over components
432 441  
433 442 // Return
434 443 return;
... ... @@ -459,52 +468,53 @@ void GModelSpectralMultiplicative::write(GXmlElement&amp; xml) const
459 468 }
460 469  
461 470 // Loop over model components
462   - for(int i = 0; i < m_spectral.size(); i++) {
  471 + for (int i = 0; i < m_spectral.size(); i++) {
463 472  
464   - // Write spectral model
465   - if (m_spectral[i] != NULL) {
  473 + // Write spectral model
  474 + if (m_spectral[i] != NULL) {
466 475  
467   - // Create new spectrum node
468   - xml.append(GXmlElement("spectrum"));
  476 + // Create new spectrum node
  477 + xml.append(GXmlElement("spectrum"));
469 478  
470   - // Get new spectrum node
471   - GXmlElement* spec = xml.element("spectrum", xml.elements("spectrum")-1);
  479 + // Get new spectrum node
  480 + GXmlElement* spec = xml.element("spectrum", xml.elements("spectrum")-1);
472 481  
473   - // Create temporary copy of the spectral model
474   - // This is a kluge to write out the original parameters.
475   - GModelSpectral* cpy = m_spectral[i]->clone();
  482 + // Create temporary copy of the spectral model. This is a kluge to
  483 + // write out the original parameters.
  484 + GModelSpectral* cpy = m_spectral[i]->clone();
476 485  
477   - // Loop over parameters of model
478   - for (int j = 0; j < cpy->size(); ++j) {
  486 + // Loop over parameters of model
  487 + for (int j = 0; j < cpy->size(); ++j) {
479 488  
480   - // Get model parameter and name
481   - GModelPar& par = (*cpy)[j];
482   - std::string parname = par.name();
  489 + // Get model parameter and name
  490 + GModelPar& par = (*cpy)[j];
  491 + std::string parname = par.name();
483 492  
484   - // Check if name contains colon
485   - if (gammalib::contains(parname, ":")) {
  493 + // Check if name contains colon
  494 + if (gammalib::contains(parname, ":")) {
486 495  
487   - // Split at the colon
488   - std::vector<std::string> splits = gammalib::split(parname, ":");
  496 + // Split at the colon
  497 + std::vector<std::string> splits = gammalib::split(parname, ":");
489 498  
490   - // Use second part of the string to recover original parameter name
491   - par.name(splits[1]);
492   - }
  499 + // Use second part of the string to recover original
  500 + // parameter name
  501 + par.name(splits[1]);
  502 + }
493 503  
494   - } // endfor: loop over parameters
  504 + } // endfor: loop over parameters
495 505  
496   - // Write spectral component
497   - cpy->write(*spec);
  506 + // Write spectral component
  507 + cpy->write(*spec);
498 508  
499   - // Add component name if previously available
500   - if (m_components[i] != gammalib::str(i+1)) {
501   - spec->attribute("component", m_components[i]);
502   - }
  509 + // Add component name if previously available
  510 + if (m_components[i] != gammalib::str(i+1)) {
  511 + spec->attribute("component", m_components[i]);
  512 + }
503 513  
504   - // Remove temporary copy
505   - delete cpy;
  514 + // Remove temporary copy
  515 + delete cpy;
506 516  
507   - } // endif: spectral model was not NULL
  517 + } // endif: spectral model was not NULL
508 518  
509 519 } // endfor: loop over model components
510 520  
... ... @@ -514,7 +524,7 @@ void GModelSpectralMultiplicative::write(GXmlElement&amp; xml) const
514 524  
515 525  
516 526 /***********************************************************************//**
517   - * @brief Print powerlaw information
  527 + * @brief Print multiplicative spectral model information
518 528 *
519 529 * @param[in] chatter Chattiness (defaults to NORMAL).
520 530 * @return String containing model information.
... ... @@ -527,19 +537,19 @@ std::string GModelSpectralMultiplicative::print(const GChatter&amp; chatter) const
527 537 // Continue only if chatter is not silent
528 538 if (chatter != SILENT) {
529 539  
530   - // Append header
  540 + // Append header
531 541 result.append("=== GModelSpectralMultiplicative ===");
532 542  
533 543 // Append information
534   - result.append("\n"+gammalib::parformat("Number of components"));
535   - result.append(gammalib::str(components()));
536   - result.append("\n"+gammalib::parformat("Number of parameters"));
537   - result.append(gammalib::str(size()));
  544 + result.append("\n"+gammalib::parformat("Number of components"));
  545 + result.append(gammalib::str(components()));
  546 + result.append("\n"+gammalib::parformat("Number of parameters"));
  547 + result.append(gammalib::str(size()));
538 548  
539   - // Print parameter information
540   - for (int i = 0; i < size(); ++i) {
541   - result.append("\n"+m_pars[i]->print(chatter));
542   - }
  549 + // Print parameter information
  550 + for (int i = 0; i < size(); ++i) {
  551 + result.append("\n"+m_pars[i]->print(chatter));
  552 + }
543 553  
544 554 } // endif: chatter was not silent
545 555  
... ... @@ -551,12 +561,13 @@ std::string GModelSpectralMultiplicative::print(const GChatter&amp; chatter) const
551 561 /***********************************************************************//**
552 562 * @brief Append spectral component
553 563 *
554   - * @param[in] spec Spectral model component
555   - * @param[in] name Name of spectral component (can be empty)
  564 + * @param[in] spec Spectral model component.
  565 + * @param[in] name Name of spectral component (can be empty).
556 566 *
557 567 * Appends a spectral component to the Multiplicative model
558 568 ***************************************************************************/
559   -void GModelSpectralMultiplicative::append(const GModelSpectral& spec, const std::string& name)
  569 +void GModelSpectralMultiplicative::append(const GModelSpectral& spec,
  570 + const std::string& name)
560 571 {
561 572 // Append model container
562 573 m_spectral.push_back(spec.clone());
... ... @@ -564,81 +575,85 @@ void GModelSpectralMultiplicative::append(const GModelSpectral&amp; spec, const std:
564 575 // Get index of latest model
565 576 int index = m_spectral.size()-1;
566 577  
567   - // Use model index if component name is empty
568   - std::string component_name = !name.empty() ? name : gammalib::str(m_spectral.size());
  578 + // Use model index as component name if component name is empty
  579 + std::string component_name = !name.empty() ? name
  580 + : gammalib::str(m_spectral.size());
  581 +
  582 + // Add component name (for now simple number)
  583 + m_components.push_back(component_name);
569 584  
570   - // Add component name (for now simple number)
571   - m_components.push_back(component_name);
  585 + // Get number of spectral parameters from model
  586 + int npars = m_spectral[index]->size();
572 587  
573   - // Get number of spectral parameters from model
574   - int npars = m_spectral[index]->size();
  588 + // Loop over model parameters
  589 + for (int ipar = 0; ipar < npars; ++ipar) {
575 590  
576   - // Loop over model parameters
577   - for (int ipar = 0; ipar < npars; ++ipar) {
  591 + // Get model parameter
  592 + GModelPar* par = &(m_spectral[index]->operator[](ipar));
578 593  
579   - // Get model parameter
580   - GModelPar* par = &(m_spectral[index]->operator[](ipar));
  594 + // Prepend component name to parameter name
  595 + par->name(component_name+":"+par->name());
581 596  
582   - // Modify parameter name
583   - par->name(component_name+":"+par->name());
  597 + // Append model parameter with new name to internal container
  598 + m_pars.push_back(par);
584 599  
585   - // Append model parameter with new name to internal container
586   - m_pars.push_back(par);
  600 + } // endfor: loop over model parameters
587 601  
588   - } // endfor: loop over model parameters
  602 + // Return
  603 + return;
589 604 }
590 605  
591 606  
592 607 /***********************************************************************//**
593   - * @brief Returns spectral component element
  608 + * @brief Return spectral model component by index
594 609 *
595 610 * @param[in] index Index of spectral component.
596   - * @return Spectral model.
  611 + * @return Pointer to spectral model.
597 612 *
598   - * Returns a spectral component to the Multiplicative model
  613 + * Returns a component of the multiplicative spectral model by @p index.
599 614 ***************************************************************************/
600 615 const GModelSpectral* GModelSpectralMultiplicative::component(const int& index) const
601 616 {
602 617 // Check if index is in validity range
603 618 if (index >= m_spectral.size() || index < 0) {
604   - throw GException::out_of_range(G_COMPONENT_INDEX, "Component Index", index, m_spectral.size(),"");
  619 + throw GException::out_of_range(G_COMPONENT_INDEX, "Component Index",
  620 + index, m_spectral.size());
605 621 }
606 622  
607 623 // Return spectral component
608 624 return m_spectral[index];
609   -
610 625 }
611 626  
612 627  
613 628 /***********************************************************************//**
614   - * @brief Returns pointer to specific spectral component
  629 + * @brief Return spectral model component by name
615 630 *
616 631 * @param[in] name Name of spectral component.
617   - * @return Spectral model.
  632 + * @return Pointer to spectral model.
618 633 *
619   - * Returns a spectral component of the Multiplicative model
  634 + * Returns a component of the multiplicative spectral model by @p name.
620 635 ***************************************************************************/
621 636 const GModelSpectral* GModelSpectralMultiplicative::component(const std::string& name) const
622 637 {
623   - // Check if model name is found
624   - int index = -1;
625   - for(int i = 0; i < m_components.size(); ++i) {
626   - if (m_components[i] == name) {
627   - index = i;
628   - break;
629   - }
630   - }
  638 + // Check if model name is found
  639 + int index = -1;
  640 + for (int i = 0; i < m_components.size(); ++i) {
  641 + if (m_components[i] == name) {
  642 + index = i;
  643 + break;
  644 + }
  645 + }
631 646  
632   - // Check if component name was found
633   - if (index == -1) {
634   - throw GException::model_not_found(G_COMPONENT_NAME, name,"");
635   - }
  647 + // Check if component name was found
  648 + if (index == -1) {
  649 + throw GException::model_not_found(G_COMPONENT_NAME, name);
  650 + }
636 651  
637 652 // Return spectral component
638 653 return m_spectral[index];
639   -
640 654 }
641 655  
  656 +
642 657 /*==========================================================================
643 658 = =
644 659 = Private methods =
... ... @@ -661,7 +676,7 @@ void GModelSpectralMultiplicative::init_members(void)
661 676 m_mc_spectrum.clear();
662 677 m_mc_emin.clear();
663 678 m_mc_emax.clear();
664   - m_n_mc_nodes = 100;
  679 + m_mc_values.clear();
665 680  
666 681 // Return
667 682 return;
... ... @@ -671,7 +686,7 @@ void GModelSpectralMultiplicative::init_members(void)
671 686 /***********************************************************************//**
672 687 * @brief Copy class members
673 688 *
674   - * @param[in] model GModelSpectralMultiplicative members which should be copied.
  689 + * @param[in] model Multiplicative spectral model.
675 690 ***************************************************************************/
676 691 void GModelSpectralMultiplicative::copy_members(const GModelSpectralMultiplicative& model)
677 692 {
... ... @@ -681,9 +696,9 @@ void GModelSpectralMultiplicative::copy_members(const GModelSpectralMultiplicati
681 696  
682 697 // Copy MC cache
683 698 m_mc_spectrum = model.m_mc_spectrum;
684   - m_n_mc_nodes = model.m_n_mc_nodes;
685 699 m_mc_emin = model.m_mc_emin;
686 700 m_mc_emax = model.m_mc_emax;
  701 + m_mc_values = model.m_mc_values;
687 702  
688 703 // Copy pointer(s) of spectral component
689 704 m_spectral.clear();
... ... @@ -691,22 +706,22 @@ void GModelSpectralMultiplicative::copy_members(const GModelSpectralMultiplicati
691 706 m_spectral.push_back(model.m_spectral[i]->clone());
692 707 }
693 708  
694   - // Store pointers to spectral parameters
  709 + // Store pointers to spectral parameters
695 710 m_pars.clear();
696   - for(int i = 0; i < model.components(); ++i) {
  711 + for (int i = 0; i < model.components(); ++i) {
697 712  
698   - // Retrieve spectral model
699   - GModelSpectral* spec = m_spectral[i];
  713 + // Retrieve spectral model
  714 + GModelSpectral* spec = m_spectral[i];
700 715  
701   - // Loop over parameters
702   - for (int ipar = 0; ipar < spec->size(); ++ipar) {
  716 + // Loop over parameters
  717 + for (int ipar = 0; ipar < spec->size(); ++ipar) {
703 718  
704   - // Get model parameter reference
705   - GModelPar& par = spec->operator[](ipar);
  719 + // Get model parameter reference
  720 + GModelPar& par = spec->operator[](ipar);
706 721  
707   - // Append model parameter pointer to internal container
708   - m_pars.push_back(&par);
709   - }
  722 + // Append model parameter pointer to internal container
  723 + m_pars.push_back(&par);
  724 + }
710 725 }
711 726  
712 727 // Return
... ... @@ -720,14 +735,15 @@ void GModelSpectralMultiplicative::copy_members(const GModelSpectralMultiplicati
720 735 void GModelSpectralMultiplicative::free_members(void)
721 736 {
722 737 // Free memory
723   - for(int i = 0; i < m_spectral.size(); ++i) {
  738 + for (int i = 0; i < m_spectral.size(); ++i) {
724 739  
725   - // Delete component i
726   - if (m_spectral[i] != NULL) delete m_spectral[i];
  740 + // Delete component i
  741 + if (m_spectral[i] != NULL) delete m_spectral[i];
727 742  
728   - // Signal free pointer
729   - m_spectral[i] = NULL;
730   - }
  743 + // Signal free pointer
  744 + m_spectral[i] = NULL;
  745 +
  746 + }
731 747  
732 748 // Return
733 749 return;
... ... @@ -743,17 +759,17 @@ void GModelSpectralMultiplicative::free_members(void)
743 759 ***************************************************************************/
744 760 void GModelSpectralMultiplicative::add_component(const GXmlElement& spec)
745 761 {
746   - // Initialise a spectral registry object
  762 + // Initialise a spectral registry object
747 763 GModelSpectralRegistry registry;
748 764  
749   - // Read spectral model
750   - GModelSpectral* ptr = registry.alloc(spec);
  765 + // Read spectral model
  766 + GModelSpectral* ptr = registry.alloc(spec);
751 767  
752   - // Get component attribute from XML file
753   - std::string component_name = spec.attribute("component");
  768 + // Get component attribute from XML file
  769 + std::string component_name = spec.attribute("component");
754 770  
755   - // Append spectral component to container
756   - append(*ptr, component_name);
  771 + // Append spectral component to container
  772 + append(*ptr, component_name);
757 773  
758 774 // Return
759 775 return;
... ... @@ -766,27 +782,66 @@ void GModelSpectralMultiplicative::add_component(const GXmlElement&amp; spec)
766 782 * @param[in] emin Minimum photon energy.
767 783 * @param[in] emax Maximum photon energy.
768 784 *
769   - * Updates the precomputation cache for Monte Carlo simulations.
  785 + * Updates the precomputation cache for Monte Carlo simulations. In case that
  786 + * the energy boundaries have changed or at least one of the model parameters
  787 + * has changed the method computes a spectral node function which has 100
  788 + * nodes per decade containing the multiplicative spectral model values and
  789 + * stores that into a Monte Carlo cache. This cache is then used by the mc()
  790 + * method for simulations.
770 791 ***************************************************************************/
771 792 void GModelSpectralMultiplicative::update_mc_cache(const GEnergy& emin,
772   - const GEnergy& emax) const
  793 + const GEnergy& emax) const
773 794  
774 795 {
775   - // Update cache if energy range has changed
776   - if (emin != m_mc_emin || emax != m_mc_emax) {
777   - m_mc_emin = emin;
778   - m_mc_emax = emax;
779   - m_mc_spectrum.clear();
780   -
781   - // Initialise energy array
782   - GEnergies energies = GEnergies(m_n_mc_nodes, m_mc_emin, m_mc_emax, true);
783   -
784   - // Append nodes to spectral function
785   - for(int i = 0; i < m_n_mc_nodes; ++i) {
786   - m_mc_spectrum.append(energies[i], eval(energies[i]));
787   - }
  796 + // Check if one of the parameters has changed. If the dimension of the
  797 + // parameter value cache differs from the number of parameters then notify
  798 + // a change. This will clear the cache and store the parameter values
  799 + // later
  800 + bool par_changed = (m_mc_values.size() != size());
  801 + if (par_changed == false) {
  802 + for (int i = 0; i < size(); ++i) {
  803 + if (m_pars[i]->value() != m_mc_values[i]) {
  804 + par_changed = true;
  805 + break;
  806 + }
  807 + }
  808 + }
788 809  
789   - } // endif: emin and emax have changed
  810 + // Update cache if energy range or parameters have changed
  811 + if (par_changed || emin != m_mc_emin || emax != m_mc_emax) {
  812 +
  813 + // Store energy range
  814 + m_mc_emin = emin;
  815 + m_mc_emax = emax;
  816 +
  817 + // If parameters have changed then store the current parameter
  818 + // values for a comparison check for the next method call
  819 + if (par_changed) {
  820 + m_mc_values.clear();
  821 + for (int i = 0; i < size(); ++i) {
  822 + m_mc_values.push_back(m_pars[i]->value());
  823 + }
  824 + }
  825 +
  826 + // Clear spectral nodes
  827 + m_mc_spectrum.clear();
  828 +
  829 + // Compute number of nodes. We use here 100 nodes per log energy and
  830 + // assure that at least 100 nodes are used.
  831 + int nodes = int((emax.log10MeV() - emin.log10MeV()) * 100.0);
  832 + if (nodes < 100) {
  833 + nodes = 100;
  834 + }
  835 +
  836 + // Initialise energy array with fixed number of nodes
  837 + GEnergies energies = GEnergies(nodes, m_mc_emin, m_mc_emax, true);
  838 +
  839 + // Append nodes to spectral function
  840 + for (int i = 0; i < energies.size(); ++i) {
  841 + m_mc_spectrum.append(energies[i], eval(energies[i]));
  842 + }
  843 +
  844 + } // endif: emin and emax have changed
790 845  
791 846 // Return
792 847 return;
... ...
test/test_GModel.cpp
  View file @eb6dc13
... ... @@ -141,7 +141,7 @@ void TestGModel::set(void)
141 141 append(static_cast<pfunction>(&TestGModel::test_logparabola),
142 142 "Test GModelSpectralLogParabola");
143 143 append(static_cast<pfunction>(&TestGModel::test_multiplicative),
144   - "Test GModelSpectralMultiplicative");
  144 + "Test GModelSpectralMultiplicative");
145 145 append(static_cast<pfunction>(&TestGModel::test_nodes),
146 146 "Test GModelSpectralNodes");
147 147 append(static_cast<pfunction>(&TestGModel::test_filefct),
... ... @@ -2020,34 +2020,35 @@ void TestGModel::test_multiplicative(void)
2020 2020 test_value(model1.type(), "Multiplicative", "Check void model type");
2021 2021  
2022 2022 // Test value constructor
2023   - model1.append(GModelSpectralPlaw(1e-17, -2.5, GEnergy(1, "TeV")), "Plaw1");
2024   - model1.append(GModelSpectralPlaw(1.0, 0.1, GEnergy(1, "TeV")), "Plaw2");
  2023 + model1.append(GModelSpectralPlaw(1.0e-17, -2.5, GEnergy(1.0, "TeV")), "Plaw1");
  2024 + model1.append(GModelSpectralPlaw(1.0, 0.1, GEnergy(1.0, "TeV")), "Plaw2");
2025 2025  
2026   - test_value(model1.eval(GEnergy(1,"TeV")),1e-17);
  2026 + // Test eval method
  2027 + test_value(model1.eval(GEnergy(1.0, "TeV")), 1.0e-17);
2027 2028 model1["Plaw2:Index"].value(2.5);
2028   - test_value(model1.eval(GEnergy(5,"TeV")),1e-17);
  2029 + test_value(model1.eval(GEnergy(5.0, "TeV")), 1.0e-17);
2029 2030  
2030 2031 // Test XML constructor
2031   - GXml xml(m_xml_point_multiplicative);
2032   - GXmlElement* element = xml.element(0)->element(0)->element("spectrum", 0);
  2032 + GXml xml(m_xml_point_multiplicative);
  2033 + GXmlElement* element = xml.element(0)->element(0)->element("spectrum", 0);
2033 2034 GModelSpectralMultiplicative model2(*element);
2034 2035 test_value(model2.size(), 7);
2035 2036 test_value(model2.type(), "Multiplicative", "Check model type");
2036   - test_value(model2.eval(GEnergy(1,"TeV")), 1e-17 * std::exp(-1));
  2037 + test_value(model2.eval(GEnergy(1.0, "TeV")), 1.0e-17 * std::exp(-1.0));
2037 2038  
2038 2039 // Create copy of model2
2039 2040 GModelSpectralMultiplicative model3(model2);
2040 2041  
2041 2042 // Change model2
2042   - model2["1:Prefactor"].value(3e-17);
2043   - test_value(model2["1:Prefactor"].value(), 3e-17);
  2043 + model2["1:Prefactor"].value(3.0e-17);
  2044 + test_value(model2["1:Prefactor"].value(), 3.0e-17);
2044 2045  
2045   - model2["2:CutoffEnergy"].value(1e7);
2046   - test_value(model2["2:CutoffEnergy"].value(), 1e7);
  2046 + model2["2:CutoffEnergy"].value(1.0e7);
  2047 + test_value(model2["2:CutoffEnergy"].value(), 1.0e7);
2047 2048  
2048   - // Verify model3 hasnt been changed
2049   - test_value(model3["1:Prefactor"].value(), 1e-17);
2050   - test_value(model3["2:CutoffEnergy"].value(), 1e6);
  2049 + // Verify model3 hasnt been changed
  2050 + test_value(model3["1:Prefactor"].value(), 1.0e-17);
  2051 + test_value(model3["2:CutoffEnergy"].value(), 1.0e6);
2051 2052  
2052 2053 // Exit test
2053 2054 return;
... ...