GammaLib  1.7.0.dev
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Macros Pages
GCTAModelSpatialMultiplicative.cpp
Go to the documentation of this file.
1 /***************************************************************************
2  * GCTAModelSpatialMultiplicative.cpp - Multiplicative spatial model class *
3  * ----------------------------------------------------------------------- *
4  * copyright (C) 2018 by Jurgen Knodlseder *
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 GCTAModelSpatialMultiplicative.cpp
23  * @brief Multiplicative spatial model class implementation
24  * @author Juergen Knoedlseder
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 "GMath.hpp"
35 #include "GRan.hpp"
36 #include "GEnergy.hpp"
37 #include "GTime.hpp"
38 #include "GXmlElement.hpp"
39 #include "GCTAObservation.hpp"
40 #include "GCTAInstDir.hpp"
41 #include "GCTARoi.hpp"
44 
45 /* __ Constants __________________________________________________________ */
46 
47 /* __ Globals ____________________________________________________________ */
49 const GCTAModelSpatialRegistry g_cta_spatial_multi_registry(&g_cta_spatial_multi_seed);
50 
51 /* __ Method name definitions ____________________________________________ */
52 #define G_READ "GCTAModelSpatialMultiplicative::read(GXmlElement&)"
53 #define G_WRITE "GCTAModelSpatialMultiplicative::write(GXmlElement&)"
54 #define G_COMPONENT_INDEX "GCTAModelSpatialMultiplicative::component(int&)"
55 #define G_COMPONENT_NAME "GCTAModelSpatialMultiplicative::component"\
56  "(std::string&)"
57 #define G_APPEND "GCTAModelSpatialMultiplicative::append(GCTAModelSpatial&, "\
58  "std::string&)"
59 
60 /* __ Macros _____________________________________________________________ */
61 
62 /* __ Coding definitions _________________________________________________ */
63 
64 /* __ Debug definitions __________________________________________________ */
65 
66 
67 /*==========================================================================
68  = =
69  = Constructors/destructors =
70  = =
71  ==========================================================================*/
72 
73 /***********************************************************************//**
74  * @brief Void constructor
75  ***************************************************************************/
78 {
79  // Initialise members
80  init_members();
81 
82  // Return
83  return;
84 }
85 
86 
87 /***********************************************************************//**
88  * @brief XML constructor
89  *
90  * @param[in] xml XML element.
91  *
92  * Creates instance of a multiplicative spatial model by extracting
93  * information from an XML element.
94  * See GCTAModelSpatialMultiplicative::read() for more information about the
95  * expected structure of the XML element.
96  ***************************************************************************/
99 {
100  // Initialise members
101  init_members();
102 
103  // Read information from XML element
104  read(xml);
105 
106  // Return
107  return;
108 }
109 
110 
111 /***********************************************************************//**
112  * @brief Copy constructor
113  *
114  * @param[in] model Spatial model.
115  ***************************************************************************/
118 {
119  // Initialise members
120  init_members();
121 
122  // Copy members
123  copy_members(model);
124 
125  // Return
126  return;
127 }
128 
129 
130 /***********************************************************************//**
131  * @brief Destructor
132  ***************************************************************************/
134 {
135  // Free members
136  free_members();
137 
138  // Return
139  return;
140 }
141 
142 
143 /*==========================================================================
144  = =
145  = Operators =
146  = =
147  ==========================================================================*/
148 
149 /***********************************************************************//**
150  * @brief Assignment operator
151  *
152  * @param[in] model Spatial model.
153  ***************************************************************************/
155 {
156  // Execute only if object is not identical
157  if (this != &model) {
158 
159  // Copy base class members
160  this->GCTAModelSpatial::operator=(model);
161 
162  // Free members
163  free_members();
164 
165  // Initialise members
166  init_members();
167 
168  // Copy members
169  copy_members(model);
170 
171  } // endif: object was not identical
172 
173  // Return
174  return *this;
175 }
176 
177 
178 /*==========================================================================
179  = =
180  = Public methods =
181  = =
182  ==========================================================================*/
183 
184 /***********************************************************************//**
185  * @brief Clear instance
186  ***************************************************************************/
188 {
189  // Free class members (base and derived classes, derived class first)
190  free_members();
192 
193  // Initialise members
195  init_members();
196 
197  // Return
198  return;
199 }
200 
201 
202 /***********************************************************************//**
203  * @brief Clone instance
204  ***************************************************************************/
206 {
207  return new GCTAModelSpatialMultiplicative(*this);
208 }
209 
210 
211 /***********************************************************************//**
212  * @brief Evaluate function (in units of sr^-1)
213  *
214  * @param[in] dir Event direction.
215  * @param[in] energy Event energy.
216  * @param[in] time Event time.
217  * @param[in] gradients Compute gradients?
218  * @return Model value
219  *
220  * Evaluate spatial model for a given event direction.
221  ***************************************************************************/
223  const GEnergy& energy,
224  const GTime& time,
225  const bool& gradients) const
226 {
227  // Initialise result
228  double value = 0.0;
229 
230  // Check for available spatial components
231  if (m_spatial.size() > 0) {
232 
233  // Initialise vector of values
234  std::vector<double> values(m_spatial.size(), 0.0);
235 
236  // Compute values of components and function value
237  for (int i = 0; i < m_spatial.size(); ++i) {
238  values[i] = m_spatial[i]->eval(dir, energy, time, gradients);
239  if (i == 0) {
240  value = values[0];
241  }
242  else {
243  value *= values[i];
244  }
245  }
246 
247  // Modify gradients if requested
248  if (gradients) {
249 
250  // Loop over model components
251  for (int i = 0; i < m_spatial.size(); ++i) {
252 
253  // Initialise scaling factor
254  double factor = 1.0;
255 
256  // Loop over other model components and compute factor
257  for (int j = 0; j < m_spatial.size(); ++j) {
258  if (i != j) {
259  factor *= values[j];
260  }
261  }
262 
263  // Loop over model parameters
264  for (int ipar = 0; ipar < m_spatial[i]->size(); ++ipar) {
265 
266  // Get reference to model parameter
267  GModelPar& par = m_spatial[i]->operator[](ipar);
268 
269  // Scale parameter gradient
270  par.gradient(par.gradient()*factor);
271 
272  } // endfor: loop over model parameters
273 
274  } // endfor: loop over models
275 
276  } //endif: compute gradients
277 
278  } // endif: there were spatial components
279 
280  // Compile option: Check for NaN/Inf
281  #if defined(G_NAN_CHECK)
282  if (gammalib::is_notanumber(value) || gammalib::is_infinite(value)) {
283  std::cout << "*** ERROR: GCTAModelSpatialMultiplicative::eval";
284  std::cout << "(dir=" << dir.print();
285  std::cout << ", energy=" << energy;
286  std::cout << ", time=" << time << "):";
287  std::cout << " NaN/Inf encountered";
288  std::cout << " (value=" << value;
289  std::cout << ")" << std::endl;
290  }
291  #endif
292 
293  // Return value
294  return value;
295 }
296 
297 
298 /***********************************************************************//**
299  * @brief Return maximum function value for Monte Carlo simulations
300  *
301  * @param[in] obs CTA Observation.
302  * @return Maximum function value for Monte Carlo simulations.
303  ***************************************************************************/
305 {
306  // Initialise maximum value
307  double value = 1.0;
308 
309  // Compute maximum value from components
310  for (int i = 0; i < m_spatial.size(); ++i) {
311  value *= m_spatial[i]->mc_max_value(obs);
312  }
313 
314  // Return maximum value
315  return value;
316 }
317 
318 
319 /***********************************************************************//**
320  * @brief Read model from XML element
321  *
322  * @param[in] xml XML element.
323  *
324  * Reads the spatial information from an XML element.
325  ***************************************************************************/
327 {
328  // Get number of spatial components
329  int n_spatials = xml.elements("spatialModel");
330 
331  // Loop over spatial elements
332  for (int i = 0; i < n_spatials; ++i) {
333 
334  // Get spatial XML element
335  const GXmlElement* spatial = xml.element("spatialModel", i);
336 
337  // Allocate a spatial registry object
338  GCTAModelSpatialRegistry registry;
339 
340  // Read spatial model
341  GCTAModelSpatial* ptr = registry.alloc(*spatial);
342 
343  // Get component attribute from XML file
344  std::string component_name = spatial->attribute("component");
345 
346  // Append spatial component to container
347  append(*ptr, component_name);
348 
349  // Free spatial model
350  delete ptr;
351 
352  } // endfor: loop over components
353 
354  // Return
355  return;
356 }
357 
358 
359 /***********************************************************************//**
360  * @brief Write model into XML element
361  *
362  * @param[in] xml XML element.
363  *
364  * @exception GException::invalid_value
365  * Existing XML element is not of the expected type.
366  *
367  * Writes the spatial information into an XML element.
368  ***************************************************************************/
370 {
371  // Set model type
372  if (xml.attribute("type") == "") {
373  xml.attribute("type", type());
374  }
375 
376  // Verify model type
377  if (xml.attribute("type") != type()) {
378  std::string msg = "Spatial model \""+xml.attribute("type")+"\" is "
379  "not of expected type \""+type()+"\".";
381  }
382 
383  // Loop over model components
384  for (int i = 0; i < m_spatial.size(); i++) {
385 
386  // Write spatial model
387  if (m_spatial[i] != NULL) {
388 
389  // Create new spectrum node
390  xml.append(GXmlElement("spatialModel"));
391 
392  // Get new spatial node
393  GXmlElement* spatial = xml.element("spatialModel",
394  xml.elements("spatialModel")-1);
395 
396  // Create temporary copy of the spatial model. This is a kluge to
397  // write out the original parameters.
398  GCTAModelSpatial* cpy = m_spatial[i]->clone();
399 
400  // Loop over parameters of model
401  for (int j = 0; j < cpy->size(); ++j) {
402 
403  // Get model parameter and name
404  GModelPar& par = (*cpy)[j];
405  std::string parname = par.name();
406 
407  // Check if name contains colon
408  if (gammalib::contains(parname, ":")) {
409 
410  // Split at the colon
411  std::vector<std::string> splits = gammalib::split(parname, ":");
412 
413  // Use second part of the string to recover original
414  // parameter name
415  par.name(splits[1]);
416  }
417 
418  } // endfor: loop over parameters
419 
420  // Write spatial component
421  cpy->write(*spatial);
422 
423  // Add component name if previously available
424  if (m_components[i] != gammalib::str(i+1)) {
425  spatial->attribute("component", m_components[i]);
426  }
427 
428  // Remove temporary copy
429  delete cpy;
430 
431  } // endif: spatial model was not NULL
432 
433  } // endfor: loop over model components
434 
435  // Return
436  return;
437 }
438 
439 
440 /***********************************************************************//**
441  * @brief Append spatial component
442  *
443  * @param[in] spatial Spatial model component.
444  * @param[in] name Name of spatial component (can be empty).
445  *
446  * @exception GException::invalid_value
447  * Invalid component name specified
448  *
449  * Appends a spatial component to the multiplicative model
450  ***************************************************************************/
452  const std::string& name)
453 {
454  // Append model container
455  m_spatial.push_back(spatial.clone());
456 
457  // Get index of latest model
458  int index = m_spatial.size()-1;
459 
460  // Use model index as component name if component name is empty
461  std::string component_name = !name.empty() ? name
462  : gammalib::str(m_spatial.size());
463 
464  // Check if component name is unique, throw exception if not
465  if (gammalib::contains(m_components, component_name)) {
466  std::string msg = "Attempt to append component with name \""+
467  component_name+"\" to multiplicative spatial model "
468  "container, but a component with the same name exists "
469  "already. Every component in the container needs a "
470  "unique name. On default the system will increment "
471  "an integer if no component name is provided.";
473  }
474 
475  // Add component name (for now simple number)
476  m_components.push_back(component_name);
477 
478  // Get number of spatial parameters from model
479  int npars = m_spatial[index]->size();
480 
481  // Loop over model parameters
482  for (int ipar = 0; ipar < npars; ++ipar) {
483 
484  // Get model parameter
485  GModelPar* par = &(m_spatial[index]->operator[](ipar));
486 
487  // Prepend component name to parameter name
488  par->name(component_name+":"+par->name());
489 
490  // Append model parameter with new name to internal container
491  m_pars.push_back(par);
492 
493  } // endfor: loop over model parameters
494 
495  // Return
496  return;
497 }
498 
499 
500 /***********************************************************************//**
501  * @brief Return spatial model component by index
502  *
503  * @param[in] index Index of spectral component.
504  * @return Pointer to spatial model.
505  *
506  * @exception GException::out_of_range
507  * Model index out of valid range
508  *
509  * Returns a component of the multiplicative spatial model by @p index.
510  ***************************************************************************/
512 {
513  // Check if index is in validity range
514  if (index >= m_spatial.size() || index < 0) {
515  throw GException::out_of_range(G_COMPONENT_INDEX, "Component Index",
516  index, m_spatial.size());
517  }
518 
519  // Return spatial component
520  return m_spatial[index];
521 }
522 
523 
524 /***********************************************************************//**
525  * @brief Return spatial model component by name
526  *
527  * @param[in] name Name of spatial component.
528  * @return Pointer to spatial model.
529  *
530  * @exception GException::invalid_argument
531  * Spatial component not found
532  *
533  * Returns a component of the multiplicative spatial model by @p name.
534  ***************************************************************************/
535 const GCTAModelSpatial* GCTAModelSpatialMultiplicative::component(const std::string& name) const
536 {
537  // Check if model name is found
538  int index = -1;
539  for (int i = 0; i < m_components.size(); ++i) {
540  if (m_components[i] == name) {
541  index = i;
542  break;
543  }
544  }
545 
546  // Check if component name was found
547  if (index == -1) {
548  std::string msg = "Model component \""+name+"\" not found in "
549  "multiplicative spatial model.";
551  }
552 
553  // Return spatial component
554  return m_spatial[index];
555 }
556 
557 
558 /***********************************************************************//**
559  * @brief Print multiplicative spatial model information
560  *
561  * @param[in] chatter Chattiness.
562  * @return String containing model information.
563  ***************************************************************************/
564 std::string GCTAModelSpatialMultiplicative::print(const GChatter& chatter) const
565 {
566  // Initialise result string
567  std::string result;
568 
569  // Continue only if chatter is not silent
570  if (chatter != SILENT) {
571 
572  // Append header
573  result.append("=== GCTAModelSpatialMultiplicative ===");
574 
575  // Append information
576  result.append("\n"+gammalib::parformat("Number of components"));
577  result.append(gammalib::str(components()));
578  result.append("\n"+gammalib::parformat("Number of parameters"));
579  result.append(gammalib::str(size()));
580 
581  // Print parameter information
582  for (int i = 0; i < size(); ++i) {
583  result.append("\n"+m_pars[i]->print(chatter));
584  }
585 
586  } // endif: chatter was not silent
587 
588  // Return result
589  return result;
590 }
591 
592 
593 /*==========================================================================
594  = =
595  = Private methods =
596  = =
597  ==========================================================================*/
598 
599 /***********************************************************************//**
600  * @brief Initialise class members
601  ***************************************************************************/
603 {
604  // Initialise model type
605  m_type = "Multiplicative";
606 
607  // Clear spatial models
608  m_spatial.clear();
609  m_components.clear();
610 
611  // Return
612  return;
613 }
614 
615 
616 /***********************************************************************//**
617  * @brief Copy class members
618  *
619  * @param[in] model Spatial background model component.
620  ***************************************************************************/
622 {
623  // Copy members
624  m_type = model.m_type;
625  m_components = model.m_components;
626 
627  // Copy pointer(s) of spatial component
628  m_spatial.clear();
629  for (int i = 0; i < model.components(); ++i) {
630  m_spatial.push_back(model.m_spatial[i]->clone());
631  }
632 
633  // Store pointers to spatial parameters
634  m_pars.clear();
635  for (int i = 0; i < model.components(); ++i) {
636 
637  // Retrieve spatial model
638  GCTAModelSpatial* spatial = m_spatial[i];
639 
640  // Loop over parameters
641  for (int ipar = 0; ipar < spatial->size(); ++ipar) {
642 
643  // Get model parameter reference
644  GModelPar& par = spatial->operator[](ipar);
645 
646  // Append model parameter pointer to internal container
647  m_pars.push_back(&par);
648  }
649  }
650 
651  // Return
652  return;
653 }
654 
655 
656 /***********************************************************************//**
657  * @brief Delete class members
658  ***************************************************************************/
660 {
661  // Free memory
662  for (int i = 0; i < m_spatial.size(); ++i) {
663 
664  // Delete component i
665  if (m_spatial[i] != NULL) delete m_spatial[i];
666 
667  // Signal free pointer
668  m_spatial[i] = NULL;
669 
670  }
671 
672  // Return
673  return;
674 }
const GCTAModelSpatialMultiplicative g_cta_spatial_multi_seed
bool contains(const std::string &str, const std::string &substring)
Checks if a substring is in a string.
Definition: GTools.cpp:1221
void free_members(void)
Delete class members.
Abstract spatial model class.
int components(void) const
Return number of spatial components.
Energy value class definition.
const std::string & name(void) const
Return parameter name.
Random number generator class definition.
XML element node class interface definition.
virtual void clear(void)
Clear instance.
double gradient(void) const
Return parameter gradient.
Multiplicative spatial model class.
virtual void write(GXmlElement &xml) const =0
virtual GCTAModelSpatial & operator=(const GCTAModelSpatial &model)
Assignment operator.
XML element node class.
Definition: GXmlElement.hpp:47
const GCTAModelSpatial * component(const int &index) const
Return spatial model component by index.
virtual double eval(const GCTAInstDir &dir, const GEnergy &energy, const GTime &time, const bool &gradients=false) const
Evaluate function (in units of sr^-1)
void copy_members(const GCTAModelSpatialMultiplicative &model)
Copy class members.
std::vector< std::string > split(const std::string &s, const std::string &sep)
Split string.
Definition: GTools.cpp:862
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.
virtual void write(GXmlElement &xml) const
Write model into XML element.
virtual double mc_max_value(const GCTAObservation &obs) const
Return maximum function value for Monte Carlo simulations.
Interface definition for the spatial model registry class.
bool is_notanumber(const double &x)
Signal if argument is not a number.
Definition: GTools.hpp:184
std::vector< std::string > m_components
Names of components.
bool is_infinite(const double &x)
Signal if argument is infinite.
Definition: GTools.hpp:167
Model parameter class.
Definition: GModelPar.hpp:87
virtual void read(const GXmlElement &xml)
Read model from XML element.
Multiplicative spatial model class interface definition.
void init_members(void)
Initialise class members.
virtual std::string print(const GChatter &chatter=NORMAL) const
Print multiplicative spatial model information.
const GXmlAttribute * attribute(const int &index) const
Return attribute.
void free_members(void)
Delete class members.
virtual GCTAModelSpatialMultiplicative & operator=(const GCTAModelSpatialMultiplicative &model)
Assignment operator.
void init_members(void)
Initialise class members.
CTA instrument direction class interface definition.
virtual GCTAModelSpatialMultiplicative * clone(void) const
Clone instance.
int size(void) const
Return number of model parameters.
virtual ~GCTAModelSpatialMultiplicative(void)
Destructor.
GChatter
Definition: GTypemaps.hpp:33
#define G_COMPONENT_NAME
virtual std::string print(const GChatter &chatter=NORMAL) const
Print instrument direction information.
#define G_COMPONENT_INDEX
CTA region of interest class interface definition.
CTA observation class interface definition.
virtual std::string type(void) const
Return model type.
virtual GXmlElement * element(const int &index)
Return pointer to GXMLElement child.
Definition: GXmlNode.cpp:634
void append(const GCTAModelSpatial &spatial, const std::string &name="")
Append spatial component.
GCTAModelSpatialMultiplicative(void)
Void constructor.
Exception handler interface definition.
CTA instrument direction class.
Definition: GCTAInstDir.hpp:59
GCTAModelSpatial * alloc(const GXmlElement &xml) const
Allocate spatial model that is found in XML element.
virtual GCTAModelSpatial * clone(void) const =0
Clones object.
virtual GXmlNode * append(const GXmlNode &node)
Append XML child node.
Definition: GXmlNode.cpp:284
CTA observation class.
std::string parformat(const std::string &s, const int &indent=0)
Convert string in parameter format.
Definition: GTools.cpp:1022
std::vector< GModelPar * > m_pars
Parameter pointers.
std::vector< GCTAModelSpatial * > m_spatial
Container of spatial models.
Time class interface definition.
Spatial model registry class definition.
Mathematical function definitions.
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