GammaLib  1.7.0.dev
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Macros Pages
GModelSpectralExponential.cpp
Go to the documentation of this file.
1 /***************************************************************************
2  * GModelSpectralExponential.cpp - Exponential spectral model class *
3  * ----------------------------------------------------------------------- *
4  * copyright (C) 2018 by Luigi Tibaldo *
5  * ----------------------------------------------------------------------- *
6  * *
7  * This program is free software: you can redistribute it and/or modify *
8  * it under the terms of the GNU General Public License as published by *
9  * the Free Software Foundation, either version 3 of the License, or *
10  * (at your option) any later version. *
11  * *
12  * This program is distributed in the hope that it will be useful, *
13  * but WITHOUT ANY WARRANTY; without even the implied warranty of *
14  * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the *
15  * GNU General Public License for more details. *
16  * *
17  * You should have received a copy of the GNU General Public License *
18  * along with this program. If not, see <http://www.gnu.org/licenses/>. *
19  * *
20  ***************************************************************************/
21 /**
22  * @file GModelSpectralExponential.cpp
23  * @brief Exponential spectral model class implementation
24  * @author Luigi Tibaldo
25  */
26 
27 /* __ Includes ___________________________________________________________ */
28 #ifdef HAVE_CONFIG_H
29 #include <config.h>
30 #endif
31 #include <cmath>
32 #include "GException.hpp"
33 #include "GTools.hpp"
34 #include "GRan.hpp"
35 #include "GIntegral.hpp"
36 #include "GEnergies.hpp"
38 #include "GModelSpectralNodes.hpp"
40 
41 /* __ Constants __________________________________________________________ */
42 
43 /* __ Globals ____________________________________________________________ */
45 const GModelSpectralRegistry g_spectral_expo_registry(&g_spectral_expo_seed);
46 
47 /* __ Method name definitions ____________________________________________ */
48 #define G_MC "GModelSpectralExponential::mc(GEnergy&, GEnergy&, GTime&, "\
49  "GRan&)"
50 #define G_WRITE "GModelSpectralExponential::write(GXmlElement&)"
51 
52 
53 /*==========================================================================
54  = =
55  = Constructors/destructors =
56  = =
57  ==========================================================================*/
58 
59 /***********************************************************************//**
60  * @brief Void constructor
61  ***************************************************************************/
64 {
65  // Initialise private members for clean destruction
66  init_members();
67 
68  // Return
69  return;
70 }
71 
72 
73 /***********************************************************************//**
74  * @brief XML constructor
75  *
76  * @param[in] xml XML element containing spectral model information.
77  *
78  * Constructs an exponential spectral model by extracting information from
79  * an XML element.
80  ***************************************************************************/
83 {
84  // Initialise members
85  init_members();
86 
87  // Read information from XML element
88  read(xml);
89 
90  // Return
91  return;
92 }
93 
94 
95 /***********************************************************************//**
96  * @brief Model constructor
97  *
98  * @param[in] spec Spectral model
99  *
100  * Constructs exponential spectral model by setting the exponent model.
101  ***************************************************************************/
104 {
105  // Initialise members
106  init_members();
107 
108  // Set exponent
109  exponent(spec);
110 
111  // Return
112  return;
113 }
114 
115 
116 /***********************************************************************//**
117  * @brief Copy constructor
118  *
119  * @param[in] model Exponential spectral model.
120  ***************************************************************************/
122  GModelSpectral(model)
123 {
124  // Initialise members
125  init_members();
126 
127  // Copy members
128  copy_members(model);
129 
130  // Return
131  return;
132 }
133 
134 
135 /***********************************************************************//**
136  * @brief Destructor
137  ***************************************************************************/
139 {
140  // Free members
141  free_members();
142 
143  // Return
144  return;
145 }
146 
147 
148 /*==========================================================================
149  = =
150  = Operators =
151  = =
152  ==========================================================================*/
153 
154 /***********************************************************************//**
155  * @brief Assignment operator
156  *
157  * @param[in] model Exponential spectral model.
158  * @return Exponential spectral model.
159  ***************************************************************************/
161 {
162  // Execute only if object is not identical
163  if (this != &model) {
164 
165  // Copy base class members
166  this->GModelSpectral::operator=(model);
167 
168  // Free members
169  free_members();
170 
171  // Initialise members
172  init_members();
173 
174  // Copy members
175  copy_members(model);
176 
177  } // endif: object was not identical
178 
179  // Return
180  return *this;
181 }
182 
183 
184 /*==========================================================================
185  = =
186  = Public methods =
187  = =
188  ==========================================================================*/
189 
190 /***********************************************************************//**
191  * @brief Clear Exponential spectral model
192  ***************************************************************************/
194 {
195  // Free class members (base and derived classes, derived class first)
196  free_members();
198 
199  // Initialise members
201  init_members();
202 
203  // Return
204  return;
205 }
206 
207 
208 /***********************************************************************//**
209  * @brief Clone Exponential spectral model model
210  *
211  * @return Pointer to deep copy of Exponential spectral model.
212  ***************************************************************************/
214 {
215  // Clone spectral power law model
216  return new GModelSpectralExponential(*this);
217 }
218 
219 
220 /***********************************************************************//**
221  * @brief Evaluate function
222  *
223  * @param[in] srcEng True photon energy.
224  * @param[in] srcTime True photon arrival time.
225  * @param[in] gradients Compute gradients?
226  * @return Model value (ph/cm2/s/MeV).
227  *
228  * Evaluates
229  *
230  * \f[
231  * \exp{M}(\rm srcEng, srcTime)
232  * \f]
233  *
234  * where \f${M}\f$ is the exponent model component.
235  *
236  * If the @p gradients flag is true the method will also compute the partial
237  * derivatives of each parameter with respect to the parameters using
238  *
239  * \f[
240  * \frac{\delta S}{\delta P_{\rm i}}\exp{M}
241  * \f]
242  *
243  * where \f${P_{\rm i}}\f$ is the i-th parameter.
244  ***************************************************************************/
246  const GTime& srcTime,
247  const bool& gradients) const
248 {
249  // Initialise result
250  double value = 0.0;
251 
252  // Check if exponent is defined
253  if (m_exponent != NULL) {
254 
255  // Calculate exponent value
256  value = m_exponent->eval(srcEng, srcTime, gradients);
257 
258  // Calculate exponential
259  value = std::exp(value);
260 
261  // Modify gradients if requested
262  if (gradients) {
263 
264  // Loop over model parameters
265  for (int ipar = 0; ipar < m_exponent->size(); ++ipar) {
266 
267  // Get reference to model parameter
268  GModelPar& par = m_exponent->operator[](ipar);
269 
270  // Scale parameter gradient
271  par.gradient(par.gradient()*value);
272 
273  } // endfor: loop over model parameters
274 
275  } // endif: compute grdients
276 
277  } //endif compute value and gradients
278 
279  // Compile option: Check for NaN/Inf
280  #if defined(G_NAN_CHECK)
281  if (gammalib::is_notanumber(value) || gammalib::is_infinite(value)) {
282  std::cout << "*** ERROR: GModelSpectralExponential::eval";
283  std::cout << "(srcEng=" << srcEng;
284  std::cout << ", srcTime=" << srcTime << "):";
285  std::cout << " NaN/Inf encountered";
286  std::cout << " (value=" << value;
287  std::cout << ")" << std::endl;
288  }
289  #endif
290 
291  // Return
292  return value;
293 }
294 
295 
296 /***********************************************************************//**
297  * @brief Returns model photon flux between [emin, emax] (units: ph/cm2/s)
298  *
299  * @param[in] emin Minimum photon energy.
300  * @param[in] emax Maximum photon energy.
301  * @return Photon flux (ph/cm2/s).
302  *
303  * Computes the photon flux of Exponential spectral model
304  ***************************************************************************/
306  const GEnergy& emax) const
307 {
308  // Initialise flux
309  double flux = 0.0;
310 
311  // Compute only if integration range is valid and exponent are available
312  if (emin < emax && m_exponent != NULL) {
313 
314  // Initialise function to integrate
315  flux_kern kernel(m_exponent);
316 
317  // Initialise integral class with function
318  GIntegral integral(&kernel);
319 
320  // Set integration precision
321  integral.eps(1.0e-8);
322 
323  // Calculate integral between emin and emax
324  flux = integral.romberg(emin.MeV(), emax.MeV());
325 
326  } // endif: integration range was valid
327 
328  // Return flux
329  return flux;
330 
331 }
332 
333 
334 /***********************************************************************//**
335  * @brief Returns model energy flux between [emin, emax] (units: erg/cm2/s)
336  *
337  * @param[in] emin Minimum photon energy.
338  * @param[in] emax Maximum photon energy.
339  * @return Energy flux (erg/cm2/s).
340  *
341  * Computes the energy flux of Exponential spectral model
342  ***************************************************************************/
344  const GEnergy& emax) const
345 {
346  // Initialise eflux
347  double eflux = 0.0;
348 
349  // Compute only if integration range is valid and exponent are available
350  if (emin < emax && m_exponent != NULL) {
351 
352  // Initialise function to integrate
353  eflux_kern kernel(m_exponent);
354 
355  // Initialise integral class with function
356  GIntegral integral(&kernel);
357 
358  // Set integration precision
359  integral.eps(1.0e-8);
360 
361  // Calculate integral between emin and emax
362  eflux = integral.romberg(emin.MeV(), emax.MeV());
363 
364  } // endif: integration range was valid
365 
366  // Return flux
367  return eflux;
368 
369 }
370 
371 
372 /***********************************************************************//**
373  * @brief Returns Monte Carlo energy between [emin, emax]
374  *
375  * @param[in] emin Minimum photon energy.
376  * @param[in] emax Maximum photon energy.
377  * @param[in] time True photon arrival time.
378  * @param[in,out] ran Random number generator.
379  * @return Energy.
380  *
381  * @exception GException::erange_invalid
382  * Energy range is invalid (emin < emax required).
383  *
384  * Returns Monte Carlo energy by randomly drawing from a Exponential
385  * spectral model.
386  ***************************************************************************/
388  const GEnergy& emax,
389  const GTime& time,
390  GRan& ran) const
391 {
392  // Throw an exception if energy range is invalid
393  if (emin >= emax) {
394  std::string msg = "Minimum energy "+emin.print()+" is equal or larger "
395  "than maximum energy "+emax.print()+". Please specify "
396  "a minimum energy that is smaller than the maximum "
397  "energy.";
399  }
400 
401  // Throw exception if exponent is undefined
402  if (m_exponent == NULL) {
403  std::string msg = "Exponent model is undefined.";
404  throw GException::invalid_value(G_MC, msg);
405  }
406 
407  // Update MC cache
408  update_mc_cache(emin, emax);
409 
410  // Set energy
411  GEnergy energy = m_mc_spectrum.mc(emin, emax, time, ran);
412 
413  // Return energy
414  return energy;
415 }
416 
417 
418 /***********************************************************************//**
419  * @brief Read model from XML element
420  *
421  * @param[in] xml XML element.
422  *
423  * Reads the spectral information from an XML element. The XML element shall
424  * have the following format
425  *
426  * <spectrum type="Exponential">
427  * <spectrum type="Constant">
428  * <parameter name="Normalization" scale="1" value="1" min="0" max="2" free="1"/>
429  * </spectrum>
430  * </spectrum>
431  *
432  ***************************************************************************/
434 {
435  // Get exponent XML element
436  const GXmlElement* spec = xml.element("spectrum", 0);
437 
438  // Continue only if the XML element contains children
439  if (spec->elements() > 0) {
440 
441  // Allocate a spectral registry object
442  GModelSpectralRegistry registry;
443 
444  // Read spectral model
445  GModelSpectral* ptr = registry.alloc(*spec);
446 
447  // Set spectral component as exponent
448  exponent(ptr);
449 
450  // Free spectral model
451  delete ptr;
452 
453  } // endif: XML element contains children
454 
455  // Return
456  return;
457 
458 }
459 
460 
461 /***********************************************************************//**
462  * @brief Write model into XML element
463  *
464  * @param[in] xml XML element.
465  *
466  * @exception GException::model_invalid_spectral
467  * Existing XML element is not of the expected type.
468  *
469  * Writes the spectral information into an XML element. The XML element
470  * will have the following format:
471  *
472  * <spectrum type="Exponential">
473  * <spectrum type="Constant">
474  * <parameter name="Normalization" scale="1" value="1" min="0" max="2" free="1"/>
475  * </spectrum>
476  * </spectrum>
477  *
478  * If no exponential model component is defined the method writes the
479  * following XML structure
480  *
481  * <spectrum type="Exponential">
482  * </spectrum>
483  *
484  ***************************************************************************/
486 {
487  // Set model type
488  if (xml.attribute("type") == "") {
489  xml.attribute("type", type());
490  }
491 
492  // Verify model type
493  if (xml.attribute("type") != type()) {
495  "Spectral model is not of type \""+type()+"\".");
496  }
497 
498  // Create a spectrum node
499  xml.append(GXmlElement("spectrum"));
500 
501  // Get spectrum node
502  GXmlElement* spec = xml.element("spectrum", 0);
503 
504  // Write spectral component if it exists
505  if (m_exponent != NULL) {
506  m_exponent->write(*spec);
507  }
508 
509  // Return
510  return;
511 }
512 
513 
514 /***********************************************************************//**
515  * @brief Set exponent
516  *
517  * @param[in] spec Spectral model to use as exponent.
518  *
519  * Set a spectral model as exponent
520  ***************************************************************************/
522 {
523  // Set exponent
524  m_exponent = spec->clone();
525 
526  // Get number of spectral parameters from model
527  int npars = m_exponent->size();
528 
529  // Store pointers to spectral parameters
530  m_pars.clear();
531  for (int ipar = 0; ipar < npars; ++ipar) {
532 
533  // Get model parameter reference
534  GModelPar& par = m_exponent->operator[](ipar);
535 
536  // Append model parameter pointer to internal container
537  m_pars.push_back(&par);
538  }
539 
540  // Return
541  return;
542 
543 }
544 
545 /***********************************************************************//**
546  * @brief Return exponent
547  *
548  * Returns pointer to exponent spectral model
549  ***************************************************************************/
551 {
552  // Returns exponent
553  return m_exponent;
554 }
555 
556 
557 /***********************************************************************//**
558  * @brief Print Exponential spectral model information
559  *
560  * @param[in] chatter Chattiness.
561  * @return String containing model information.
562  ***************************************************************************/
563 std::string GModelSpectralExponential::print(const GChatter& chatter) const
564 {
565  // Initialise result string
566  std::string result;
567 
568  // Continue only if chatter is not silent
569  if (chatter != SILENT) {
570 
571  // Append header
572  result.append("=== GModelSpectralExponential ===");
573 
574  // Append information
575  result.append("\n"+gammalib::parformat("Number of parameters"));
576  result.append(gammalib::str(size()));
577 
578  // Print parameter information
579  for (int i = 0; i < size(); ++i) {
580  result.append("\n"+m_pars[i]->print(chatter));
581  }
582 
583  } // endif: chatter was not silent
584 
585  // Return result
586  return result;
587 }
588 
589 
590 /*==========================================================================
591  = =
592  = Private methods =
593  = =
594  ==========================================================================*/
595 
596 /***********************************************************************//**
597  * @brief Initialise class members
598  ***************************************************************************/
600 {
601  // Initialise model type
602  m_type = "Exponential";
603 
604  // Clear exponent
605  m_exponent = NULL;
606 
607  // Clear MC cache
609  m_mc_emin.clear();
610  m_mc_emax.clear();
611  m_mc_values.clear();
612 
613  // Return
614  return;
615 }
616 
617 
618 /***********************************************************************//**
619  * @brief Copy class members
620  *
621  * @param[in] model Exponential spectral model.
622  ***************************************************************************/
624 {
625  // Copy members
626  m_type = model.m_type;
627 
628  // Copy MC cache
630  m_mc_emin = model.m_mc_emin;
631  m_mc_emax = model.m_mc_emax;
632  m_mc_values = model.m_mc_values;
633 
634  // Set exponent
635  if (model.m_exponent != NULL) {
636  exponent(model.m_exponent);
637  }
638 
639  // Return
640  return;
641 }
642 
643 
644 /***********************************************************************//**
645  * @brief Delete class members
646  ***************************************************************************/
648 {
649  // Free exponent
650  if (m_exponent != NULL) delete m_exponent;
651 
652  // Signal free pointer
653  m_exponent = NULL;
654 
655  // Return
656  return;
657 }
658 
659 
660 /***********************************************************************//**
661  * @brief Update Monte Carlo pre computation cache
662  *
663  * @param[in] emin Minimum photon energy.
664  * @param[in] emax Maximum photon energy.
665  *
666  * Updates the precomputation cache for Monte Carlo simulations. In case that
667  * the energy boundaries have changed or at least one of the model parameters
668  * has changed the method computes a spectral node function which has 100
669  * nodes per decade containing the Exponential spectral model values and
670  * stores that into a Monte Carlo cache. This cache is then used by the mc()
671  * method for simulations.
672  ***************************************************************************/
674  const GEnergy& emax) const
675 
676 {
677  // Check if one of the parameters has changed. If the dimension of the
678  // parameter value cache differs from the number of parameters then notify
679  // a change. This will clear the cache and store the parameter values
680  // later
681  bool par_changed = (m_mc_values.size() != size());
682  if (par_changed == false) {
683  for (int i = 0; i < size(); ++i) {
684  if (m_pars[i]->value() != m_mc_values[i]) {
685  par_changed = true;
686  break;
687  }
688  }
689  }
690 
691  // Update cache if energy range or parameters have changed
692  if (par_changed || emin != m_mc_emin || emax != m_mc_emax) {
693 
694  // Store energy range
695  m_mc_emin = emin;
696  m_mc_emax = emax;
697 
698  // If parameters have changed then store the current parameter
699  // values for a comparison check for the next method call
700  if (par_changed) {
701  m_mc_values.clear();
702  for (int i = 0; i < size(); ++i) {
703  m_mc_values.push_back(m_pars[i]->value());
704  }
705  }
706 
707  // Clear spectral nodes
709 
710  // Compute number of nodes. We use here 100 nodes per log energy and
711  // assure that at least 100 nodes are used.
712  int nodes = int((emax.log10MeV() - emin.log10MeV()) * 100.0);
713  if (nodes < 100) {
714  nodes = 100;
715  }
716 
717  // Initialise energy array with fixed number of nodes
718  GEnergies energies = GEnergies(nodes, m_mc_emin, m_mc_emax, true);
719 
720  // Append nodes to spectral function
721  for (int i = 0; i < energies.size(); ++i) {
722  m_mc_spectrum.append(energies[i], eval(energies[i]));
723  }
724 
725  } // endif: emin and emax have changed
726 
727  // Return
728  return;
729 }
virtual void write(GXmlElement &xml) const =0
virtual void write(GXmlElement &xml) const
Write model into XML element.
void free_members(void)
Delete class members.
void update_mc_cache(const GEnergy &emin, const GEnergy &emax) const
Update Monte Carlo pre computation cache.
#define G_WRITE
virtual GEnergy mc(const GEnergy &emin, const GEnergy &emax, const GTime &time, GRan &ran) const
Returns MC energy between [emin, emax].
double romberg(std::vector< double > bounds, const int &order=5)
Perform Romberg integration.
Definition: GIntegral.cpp:380
Random number generator class definition.
Abstract spectral model base class.
GModelSpectral * m_exponent
Exponent model component.
void init_members(void)
Initialise class members.
double gradient(void) const
Return parameter gradient.
int size(void) const
Return number of parameters.
virtual GModelSpectral & operator=(const GModelSpectral &model)
Assignment operator.
virtual std::string type(void) const
Return model type.
GEnergy m_mc_emin
Last minimum energy.
virtual double eval(const GEnergy &srcEng, const GTime &srcTime=GTime(), const bool &gradients=false) const =0
virtual GModelSpectralExponential * clone(void) const
Clone Exponential spectral model model.
virtual void read(const GXmlElement &xml)
Read model from XML element.
std::vector< GModelPar * > m_pars
Parameter pointers.
XML element node class.
Definition: GXmlElement.hpp:47
Spectral model registry class definition.
Exponential spectral model class interface definition.
Random number generator class.
Definition: GRan.hpp:44
Spectral nodes model class definition.
double MeV(void) const
Return energy in MeV.
Definition: GEnergy.cpp:321
virtual int elements(void) const
Return number of GXMLElement children of node.
Definition: GXmlNode.cpp:580
Time class.
Definition: GTime.hpp:54
Gammalib tools definition.
GIntegral class interface definition.
Definition: GIntegral.hpp:45
const GModelSpectralExponential g_spectral_expo_seed
virtual void clear(void)
Clear Exponential spectral model.
Energy container class.
Definition: GEnergies.hpp:60
double log10MeV(void) const
Return log10 of energy in MeV.
Definition: GEnergy.cpp:423
virtual double eval(const GEnergy &srcEng, const GTime &srcTime=GTime(), const bool &gradients=false) const
Evaluate function.
bool is_notanumber(const double &x)
Signal if argument is not a number.
Definition: GTools.hpp:184
bool is_infinite(const double &x)
Signal if argument is infinite.
Definition: GTools.hpp:167
void append(const GEnergy &energy, const double &intensity)
Append node.
Model parameter class.
Definition: GModelPar.hpp:87
const GXmlAttribute * attribute(const int &index) const
Return attribute.
std::string print(const GChatter &chatter=NORMAL) const
Print energy.
Definition: GEnergy.cpp:748
Energy container class definition.
int size(void) const
Return number of energies in container.
Definition: GEnergies.hpp:162
virtual GModelSpectralExponential & operator=(const GModelSpectralExponential &model)
Assignment operator.
virtual double flux(const GEnergy &emin, const GEnergy &emax) const
Returns model photon flux between emin, emax
GModelSpectral * alloc(const GXmlElement &xml) const
Allocate spectral model that is found in XML element.
GModelSpectralNodes m_mc_spectrum
MC spectrum cache.
const GModelSpectral * exponent(void) const
Return exponent.
GChatter
Definition: GTypemaps.hpp:33
Exponential spectral model class.
#define G_MC
Interface definition for the spectral model registry class.
virtual GModelSpectral * clone(void) const =0
Clones object.
void eps(const double &eps)
Set relative precision.
Definition: GIntegral.hpp:206
virtual GXmlElement * element(const int &index)
Return pointer to GXMLElement child.
Definition: GXmlNode.cpp:634
virtual void clear(void)
Clear spectral nodes model.
void init_members(void)
Initialise class members.
virtual ~GModelSpectralExponential(void)
Destructor.
Exception handler interface definition.
GVector exp(const GVector &vector)
Computes exponential of vector elements.
Definition: GVector.cpp:1142
GEnergy m_mc_emax
Last maximum energy.
virtual double eflux(const GEnergy &emin, const GEnergy &emax) const
Returns model energy flux between emin, emax
virtual GXmlNode * append(const GXmlNode &node)
Append XML child node.
Definition: GXmlNode.cpp:284
std::string parformat(const std::string &s, const int &indent=0)
Convert string in parameter format.
Definition: GTools.cpp:1022
virtual GEnergy mc(const GEnergy &emin, const GEnergy &emax, const GTime &time, GRan &ran) const
Returns Monte Carlo energy between [emin, emax].
Integration class interface definition.
GModelSpectralExponential(void)
Void constructor.
void copy_members(const GModelSpectralExponential &model)
Copy class members.
void free_members(void)
Delete class members.
void clear(void)
Clear instance.
Definition: GEnergy.cpp:261
std::vector< double > m_mc_values
Parameter values.
Class that handles energies in a unit independent way.
Definition: GEnergy.hpp:48
std::string str(const unsigned short int &value)
Convert unsigned short integer value into string.
Definition: GTools.cpp:413
virtual std::string print(const GChatter &chatter=NORMAL) const
Print Exponential spectral model information.