GammaLib  2.1.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-2021 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  * Returns Monte Carlo energy by randomly drawing from a Exponential
382  * spectral model.
383  ***************************************************************************/
385  const GEnergy& emax,
386  const GTime& time,
387  GRan& ran) const
388 {
389  // Check energy interval
391 
392  // Throw exception if exponent is undefined
393  if (m_exponent == NULL) {
394  std::string msg = "Exponent model is undefined.";
395  throw GException::invalid_value(G_MC, msg);
396  }
397 
398  // Update MC cache
399  update_mc_cache(emin, emax);
400 
401  // Set energy
402  GEnergy energy = m_mc_spectrum.mc(emin, emax, time, ran);
403 
404  // Return energy
405  return energy;
406 }
407 
408 
409 /***********************************************************************//**
410  * @brief Read model from XML element
411  *
412  * @param[in] xml XML element.
413  *
414  * Reads the spectral information from an XML element. The XML element shall
415  * have the following format
416  *
417  * <spectrum type="Exponential">
418  * <spectrum type="Constant">
419  * <parameter name="Normalization" scale="1" value="1" min="0" max="2" free="1"/>
420  * </spectrum>
421  * </spectrum>
422  *
423  ***************************************************************************/
425 {
426  // Get exponent XML element
427  const GXmlElement* spec = xml.element("spectrum", 0);
428 
429  // Continue only if the XML element contains children
430  if (spec->elements() > 0) {
431 
432  // Allocate a spectral registry object
433  GModelSpectralRegistry registry;
434 
435  // Read spectral model
436  GModelSpectral* ptr = registry.alloc(*spec);
437 
438  // Set spectral component as exponent
439  exponent(ptr);
440 
441  // Free spectral model
442  delete ptr;
443 
444  } // endif: XML element contains children
445 
446  // Return
447  return;
448 
449 }
450 
451 
452 /***********************************************************************//**
453  * @brief Write model into XML element
454  *
455  * @param[in] xml XML element.
456  *
457  * Writes the spectral information into an XML element. The XML element
458  * will have the following format:
459  *
460  * <spectrum type="Exponential">
461  * <spectrum type="Constant">
462  * <parameter name="Normalization" scale="1" value="1" min="0" max="2" free="1"/>
463  * </spectrum>
464  * </spectrum>
465  *
466  * If no exponential model component is defined the method writes the
467  * following XML structure
468  *
469  * <spectrum type="Exponential">
470  * </spectrum>
471  *
472  ***************************************************************************/
474 {
475  // Verify model type
477 
478  // Create a spectrum node
479  xml.append(GXmlElement("spectrum"));
480 
481  // Get spectrum node
482  GXmlElement* spec = xml.element("spectrum", 0);
483 
484  // Write spectral component if it exists
485  if (m_exponent != NULL) {
486  m_exponent->write(*spec);
487  }
488 
489  // Return
490  return;
491 }
492 
493 
494 /***********************************************************************//**
495  * @brief Set exponent
496  *
497  * @param[in] spec Spectral model to use as exponent.
498  *
499  * Set a spectral model as exponent
500  ***************************************************************************/
502 {
503  // Set exponent
504  m_exponent = spec->clone();
505 
506  // Get number of spectral parameters from model
507  int npars = m_exponent->size();
508 
509  // Store pointers to spectral parameters
510  m_pars.clear();
511  for (int ipar = 0; ipar < npars; ++ipar) {
512 
513  // Get model parameter reference
514  GModelPar& par = m_exponent->operator[](ipar);
515 
516  // Append model parameter pointer to internal container
517  m_pars.push_back(&par);
518  }
519 
520  // Return
521  return;
522 
523 }
524 
525 /***********************************************************************//**
526  * @brief Return exponent
527  *
528  * Returns pointer to exponent spectral model
529  ***************************************************************************/
531 {
532  // Returns exponent
533  return m_exponent;
534 }
535 
536 
537 /***********************************************************************//**
538  * @brief Print Exponential spectral model information
539  *
540  * @param[in] chatter Chattiness.
541  * @return String containing model information.
542  ***************************************************************************/
543 std::string GModelSpectralExponential::print(const GChatter& chatter) const
544 {
545  // Initialise result string
546  std::string result;
547 
548  // Continue only if chatter is not silent
549  if (chatter != SILENT) {
550 
551  // Append header
552  result.append("=== GModelSpectralExponential ===");
553 
554  // Append information
555  result.append("\n"+gammalib::parformat("Number of parameters"));
556  result.append(gammalib::str(size()));
557 
558  // Print parameter information
559  for (int i = 0; i < size(); ++i) {
560  result.append("\n"+m_pars[i]->print(chatter));
561  }
562 
563  } // endif: chatter was not silent
564 
565  // Return result
566  return result;
567 }
568 
569 
570 /*==========================================================================
571  = =
572  = Private methods =
573  = =
574  ==========================================================================*/
575 
576 /***********************************************************************//**
577  * @brief Initialise class members
578  ***************************************************************************/
580 {
581  // Initialise model type
582  m_type = "Exponential";
583 
584  // Clear exponent
585  m_exponent = NULL;
586 
587  // Clear MC cache
589  m_mc_emin.clear();
590  m_mc_emax.clear();
591  m_mc_values.clear();
592 
593  // Return
594  return;
595 }
596 
597 
598 /***********************************************************************//**
599  * @brief Copy class members
600  *
601  * @param[in] model Exponential spectral model.
602  ***************************************************************************/
604 {
605  // Copy members
606  m_type = model.m_type;
607 
608  // Copy MC cache
610  m_mc_emin = model.m_mc_emin;
611  m_mc_emax = model.m_mc_emax;
612  m_mc_values = model.m_mc_values;
613 
614  // Set exponent
615  if (model.m_exponent != NULL) {
616  exponent(model.m_exponent);
617  }
618 
619  // Return
620  return;
621 }
622 
623 
624 /***********************************************************************//**
625  * @brief Delete class members
626  ***************************************************************************/
628 {
629  // Free exponent
630  if (m_exponent != NULL) delete m_exponent;
631 
632  // Signal free pointer
633  m_exponent = NULL;
634 
635  // Return
636  return;
637 }
638 
639 
640 /***********************************************************************//**
641  * @brief Update Monte Carlo pre computation cache
642  *
643  * @param[in] emin Minimum photon energy.
644  * @param[in] emax Maximum photon energy.
645  *
646  * Updates the precomputation cache for Monte Carlo simulations. In case that
647  * the energy boundaries have changed or at least one of the model parameters
648  * has changed the method computes a spectral node function which has 100
649  * nodes per decade containing the Exponential spectral model values and
650  * stores that into a Monte Carlo cache. This cache is then used by the mc()
651  * method for simulations.
652  ***************************************************************************/
654  const GEnergy& emax) const
655 
656 {
657  // Check if one of the parameters has changed. If the dimension of the
658  // parameter value cache differs from the number of parameters then notify
659  // a change. This will clear the cache and store the parameter values
660  // later
661  bool par_changed = (m_mc_values.size() != size());
662  if (par_changed == false) {
663  for (int i = 0; i < size(); ++i) {
664  if (m_pars[i]->value() != m_mc_values[i]) {
665  par_changed = true;
666  break;
667  }
668  }
669  }
670 
671  // Update cache if energy range or parameters have changed
672  if (par_changed || emin != m_mc_emin || emax != m_mc_emax) {
673 
674  // Store energy range
675  m_mc_emin = emin;
676  m_mc_emax = emax;
677 
678  // If parameters have changed then store the current parameter
679  // values for a comparison check for the next method call
680  if (par_changed) {
681  m_mc_values.clear();
682  for (int i = 0; i < size(); ++i) {
683  m_mc_values.push_back(m_pars[i]->value());
684  }
685  }
686 
687  // Clear spectral nodes
689 
690  // Compute number of nodes. We use here 100 nodes per log energy and
691  // assure that at least 100 nodes are used.
692  int nodes = int((emax.log10MeV() - emin.log10MeV()) * 100.0);
693  if (nodes < 100) {
694  nodes = 100;
695  }
696 
697  // Initialise energy array with fixed number of nodes
698  GEnergies energies = GEnergies(nodes, m_mc_emin, m_mc_emax, "LOG");
699 
700  // Append nodes to spectral function
701  for (int i = 0; i < energies.size(); ++i) {
702  m_mc_spectrum.append(energies[i], eval(energies[i]));
703  }
704 
705  } // endif: emin and emax have changed
706 
707  // Return
708  return;
709 }
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:382
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:48
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:586
Time class.
Definition: GTime.hpp:55
Gammalib tools definition.
GIntegral class interface definition.
Definition: GIntegral.hpp:46
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:201
bool is_infinite(const double &x)
Signal if argument is infinite.
Definition: GTools.hpp:184
void append(const GEnergy &energy, const double &intensity)
Append node.
Model parameter class.
Definition: GModelPar.hpp:87
Energy container class definition.
int size(void) const
Return number of energies in container.
Definition: GEnergies.hpp:170
virtual GModelSpectralExponential & operator=(const GModelSpectralExponential &model)
Assignment operator.
void check_energy_interval(const std::string &origin, const GEnergy &emin, const GEnergy &emax)
Checks energy interval.
Definition: GException.cpp:364
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:214
virtual GXmlElement * element(const int &index)
Return pointer to GXMLElement child.
Definition: GXmlNode.cpp:640
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:1232
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:287
std::string parformat(const std::string &s, const int &indent=0)
Convert string in parameter format.
Definition: GTools.cpp:1143
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:489
virtual std::string print(const GChatter &chatter=NORMAL) const
Print Exponential spectral model information.
void xml_check_type(const std::string &origin, GXmlElement &xml, const std::string &type)
Checks the model type.
Definition: GTools.cpp:1819