GammaLib  2.0.0
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Macros Pages
GObservations.cpp
Go to the documentation of this file.
1 /***************************************************************************
2  * GObservations.cpp - Observation container class *
3  * ----------------------------------------------------------------------- *
4  * copyright (C) 2009-2021 by Juergen Knoedlseder *
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 GObservations.cpp
23  * @brief Observations container class implementation
24  * @author Juergen Knoedlseder
25  */
26 
27 /* __ Includes ___________________________________________________________ */
28 #ifdef HAVE_CONFIG_H
29 #include <config.h>
30 #endif
31 #include "GTools.hpp"
32 #include "GException.hpp"
33 #include "GFilename.hpp"
34 #include "GXml.hpp"
35 #include "GModel.hpp"
36 #include "GObservations.hpp"
37 #include "GObservationRegistry.hpp"
38 #include "GMatrixSparse.hpp"
39 #include "GOptimizer.hpp"
40 
41 /* __ Method name definitions ____________________________________________ */
42 #define G_AT "GObservations::at(int&)"
43 #define G_SET "GObservations::set(int&, GObservation&)"
44 #define G_APPEND "GObservations::append(GObservation&)"
45 #define G_INSERT "GObservations::insert(int&, GObservation&)"
46 #define G_REMOVE "GObservations::remove(int&)"
47 #define G_EXTEND "GObservations::extend(GObservations&)"
48 #define G_READ "GObservations::read(GXml&)"
49 #define G_NPRED "GObservations::npred(std::string&, std::string&)"
50 
51 /* __ Macros _____________________________________________________________ */
52 
53 /* __ Coding definitions _________________________________________________ */
54 
55 /* __ Debug definitions __________________________________________________ */
56 
57 
58 /*==========================================================================
59  = =
60  = Constructors/destructors =
61  = =
62  ==========================================================================*/
63 
64 /***********************************************************************//**
65  * @brief Void constructor
66  ***************************************************************************/
68 {
69  // Initialise members
70  init_members();
71 
72  // Return
73  return;
74 }
75 
76 
77 /***********************************************************************//**
78  * @brief Copy constructor
79  *
80  * @param obs Observation container.
81  ***************************************************************************/
83 {
84  // Initialise members
85  init_members();
86 
87  // Copy members
88  copy_members(obs);
89 
90  // Return
91  return;
92 }
93 
94 
95 /***********************************************************************//**
96  * @brief Load constructor
97  *
98  * @param[in] filename XML filename.
99  *
100  * Construct the observation container by loading all relevant information
101  * from a XML file. Please refer to the read() method for more information
102  * about the structure of the XML file.
103  ***************************************************************************/
105 {
106  // Initialise members
107  init_members();
108 
109  // Load XML file
110  load(filename);
111 
112  // Return
113  return;
114 }
115 
116 
117 /***********************************************************************//**
118  * @brief Destructor
119  ***************************************************************************/
121 {
122  // Free members
123  free_members();
124 
125  // Return
126  return;
127 }
128 
129 
130 /*==========================================================================
131  = =
132  = Operators =
133  = =
134  ==========================================================================*/
135 
136 /***********************************************************************//**
137  * @brief Assignment operator
138  *
139  * @param[in] obs Observation container.
140  * @return Observation container.
141  ***************************************************************************/
143 {
144  // Execute only if object is not identical
145  if (this != &obs) {
146 
147  // Free members
148  free_members();
149 
150  // Initialise members
151  init_members();
152 
153  // Copy members
154  copy_members(obs);
155 
156  } // endif: object was not identical
157 
158  // Return this object
159  return *this;
160 }
161 
162 
163 /*==========================================================================
164  = =
165  = Public methods =
166  = =
167  ==========================================================================*/
168 
169 /***********************************************************************//**
170  * @brief Clear observations
171  ***************************************************************************/
173 {
174  // Free members
175  free_members();
176 
177  // Initialise members
178  init_members();
179 
180  // Return
181  return;
182 }
183 
184 
185 /***********************************************************************//**
186  * @brief Clone observations
187  *
188  * @return Pointer to deep copy of observation container.
189  ***************************************************************************/
191 {
192  // Clone observations
193  return new GObservations(*this);
194 }
195 
196 
197 /***********************************************************************//**
198  * @brief Return pointer to observation
199  *
200  * @param[in] index Observation index [0,...,size()[.
201  * @return Pointer to observation.
202  *
203  * @exception GException::out_of_range
204  * Operation index is out of range.
205  *
206  * Returns a pointer to the observation with specified @p index.
207  ***************************************************************************/
208 GObservation* GObservations::at(const int& index)
209 {
210  // Raise exception if index is out of range
211  if (index < 0 || index >= size()) {
212  throw GException::out_of_range(G_AT, "Observation index",
213  index, size());
214  }
215 
216  // Return pointer
217  return m_obs[index];
218 }
219 
220 
221 /***********************************************************************//**
222  * @brief Return pointer to observation (const version)
223  *
224  * @param[in] index Observation index [0,...,size()[.
225  * @return Pointer to observation.
226  *
227  * @exception GException::out_of_range
228  * Operation index is out of range.
229  *
230  * Returns a const pointer to the observation with specified @p index.
231  ***************************************************************************/
232 const GObservation* GObservations::at(const int& index) const
233 {
234  // Raise exception if index is out of range
235  if (index < 0 || index >= size()) {
236  throw GException::out_of_range(G_AT, "Observation index",
237  index, size());
238  }
239 
240  // Return pointer
241  return m_obs[index];
242 }
243 
244 
245 /***********************************************************************//**
246  * @brief Set observation in container
247  *
248  * @param[in] index Observation index [0,...,size()[.
249  * @param[in] obs Observation.
250  * @return Pointer to deep copy of observation.
251  *
252  * @exception GException::out_of_range
253  * Observation index is out of range.
254  * @exception GException::invalid_value
255  * Observation with same instrument and identifier already
256  * exists in container.
257  *
258  * Set a deep copy and observation @p obs at the specified @p index in the
259  * container.
260  ***************************************************************************/
261 GObservation* GObservations::set(const int& index, const GObservation& obs)
262 {
263  // Compile option: raise exception if index is out of range
264  #if defined(G_RANGE_CHECK)
265  if (index < 0 || index >= size()) {
266  throw GException::out_of_range(G_SET, "Observation index",
267  index, size());
268  }
269  #endif
270 
271  // Raise an exception if an observation with specified instrument and
272  // identifier already exists
273  int inx = get_index(obs.instrument(), obs.id());
274  if (inx != -1 && inx != index) {
275  std::string msg =
276  "Attempt to set \""+obs.instrument()+"\" observation with"
277  " identifier \""+obs.id()+"\" in observation container at"
278  " index "+gammalib::str(index)+", but an observation with the"
279  " same attributes exists already at index "+gammalib::str(inx)+
280  " in the container.\n"
281  "Every observation for a given instrument in the observation"
282  " container needs a unique identifier.";
283  throw GException::invalid_value(G_SET, msg);
284  }
285 
286  // Delete existing observation
287  if (m_obs[index] != NULL) delete m_obs[index];
288 
289  // Store pointer to a deep copy of the observation
290  m_obs[index] = obs.clone();
291 
292  // Return pointer
293  return m_obs[index];
294 }
295 
296 
297 /***********************************************************************//**
298  * @brief Append observation to container
299  *
300  * @param[in] obs Observation.
301  * @return Pointer to deep copy of observation.
302  *
303  * @exception GException::invalid_value
304  * Observation with same instrument and identifier already
305  * exists in container.
306  *
307  * Appends a deep copy of an observation to the container.
308  ***************************************************************************/
310 {
311  // Raise an exception if an observation with specified instrument and
312  // identifier already exists
313  int inx = get_index(obs.instrument(), obs.id());
314  if (inx != -1) {
315  std::string msg =
316  "Attempt to append \""+obs.instrument()+"\" observation with"
317  " identifier \""+obs.id()+"\" to observation container, but an"
318  " observation with the same attributes exists already at"
319  " index "+gammalib::str(inx)+" in the container.\n"
320  "Every observation for a given instrument in the observation"
321  " container needs a unique identifier.";
323  }
324 
325  // Clone observation
326  GObservation* ptr = obs.clone();
327 
328  // Append to list
329  m_obs.push_back(ptr);
330 
331  // Return pointer
332  return ptr;
333 }
334 
335 
336 /***********************************************************************//**
337  * @brief Insert observation into container
338  *
339  * @param[in] index Observation index [0,...,size()[.
340  * @param[in] obs Observation.
341  * @return Pointer to deep copy of observation.
342  *
343  * @exception GException::out_of_range
344  * Observation index is out of range.
345  * @exception GException::invalid_value
346  * Observation with same instrument and identifier already
347  * exists in container.
348  *
349  * Inserts a deep copy of an observation into the container before the
350  * observation with the specified @p index.
351  ***************************************************************************/
352 GObservation* GObservations::insert(const int& index, const GObservation& obs)
353 {
354  // Compile option: raise exception if index is out of range
355  #if defined(G_RANGE_CHECK)
356  if (is_empty()) {
357  if (index > 0) {
358  throw GException::out_of_range(G_INSERT, "Observation index",
359  index, size());
360  }
361  }
362  else {
363  if (index < 0 || index >= size()) {
364  throw GException::out_of_range(G_INSERT, "Observation index",
365  index, size());
366  }
367  }
368  #endif
369 
370  // Raise an exception if an observation with specified instrument and
371  // identifier already exists
372  int inx = get_index(obs.instrument(), obs.id());
373  if (inx != -1 && inx != index) {
374  std::string msg =
375  "Attempt to insert \""+obs.instrument()+"\" observation with"
376  " identifier \""+obs.id()+"\" in observation container before"
377  " index "+gammalib::str(index)+", but an observation with the"
378  " same attributes exists already at index "+gammalib::str(inx)+
379  " in the container.\n"
380  "Every observation for a given instrument in the observation"
381  " container needs a unique identifier.";
383  }
384 
385  // Clone observation
386  GObservation* ptr = obs.clone();
387 
388  // Clone observation and insert into list
389  m_obs.insert(m_obs.begin()+index, ptr);
390 
391  // Return pointer
392  return ptr;
393 }
394 
395 
396 /***********************************************************************//**
397  * @brief Remove observation from container
398  *
399  * @param[in] index Observation index [0,...,size()[.
400  *
401  * @exception GException::out_of_range
402  * Observation index is out of range.
403  *
404  * Removes observation of specified @p index from the container.
405  ***************************************************************************/
406 void GObservations::remove(const int& index)
407 {
408  // Compile option: If index is outside boundary then raise exception
409  #if defined(G_RANGE_CHECK)
410  if (index < 0 || index >= size()) {
411  throw GException::out_of_range(G_REMOVE, "Observation index",
412  index, size());
413  }
414  #endif
415 
416  // Delete observation
417  delete m_obs[index];
418 
419  // Erase observation component from container
420  m_obs.erase(m_obs.begin() + index);
421 
422  // Return
423  return;
424 }
425 
426 
427 /***********************************************************************//**
428  * @brief Append observations from observation container
429  *
430  * @param[in] obs Observations.
431  *
432  * @exception GException::invalid_value
433  * Observation with same instrument and identifier already
434  * exists in container.
435  *
436  * Appends deep copies of observations to the container.
437  ***************************************************************************/
439 {
440  // Do nothing if observation container is empty
441  if (!obs.is_empty()) {
442 
443  // Get size. Note that we extract the size first to avoid an
444  // endless loop that arises when a container is appended to
445  // itself.
446  int num = obs.size();
447 
448  // Reserve enough space
449  reserve(size() + num);
450 
451  // Loop over all observations, clone them and append them to the
452  // list
453  for (int i = 0; i < num; ++i) {
454 
455  // Raise an exception if an observation with specified
456  // instrument and identifier already exists
457  int inx = get_index(obs[i]->instrument(), obs[i]->id());
458  if (inx != -1) {
459  std::string msg =
460  "Attempt to append \""+obs[i]->instrument()+"\""
461  " observation with identifier \""+obs[i]->id()+"\""
462  " to observation container, but an observation with the"
463  " same attributes exists already at index "+
464  gammalib::str(inx)+" in the container.\n"
465  "Every observation for a given instrument in the"
466  " observation container needs a unique identifier.";
468  }
469 
470  // Append observation to container
471  m_obs.push_back(obs[i]->clone());
472  }
473 
474  } // endif: observation container was not empty
475 
476  // Return
477  return;
478 }
479 
480 
481 /***********************************************************************//**
482  * @brief Signals if observation exists
483  *
484  * @param[in] instrument Instrument.
485  * @param[in] id Observation identifier.
486  * @return True if observation with specified @p instrument and identifier
487  * @p id exists.
488  *
489  * Searches all observations for a match with the specified @p instrument
490  * and @p identifier. If the specified attributes have been found, true is
491  * returned.
492  ***************************************************************************/
493 bool GObservations::contains(const std::string& instrument,
494  const std::string& id) const
495 {
496  // Get observation index
497  int index = get_index(instrument, id);
498 
499  // Return
500  return (index != -1);
501 }
502 
503 
504 /***********************************************************************//**
505  * @brief Load observations from XML file
506  *
507  * @param[in] filename XML filename.
508  *
509  * Loads observation from a XML file into the container. Please refer to the
510  * read() method for more information about the structure of the XML file.
511  ***************************************************************************/
512 void GObservations::load(const GFilename& filename)
513 {
514  // Clear any existing observations
515  clear();
516 
517  // Load XML document
518  GXml xml(filename.url());
519 
520  // Read observations from XML document
521  read(xml);
522 
523  // Return
524  return;
525 }
526 
527 
528 /***********************************************************************//**
529  * @brief Save observations into XML file
530  *
531  * @param[in] filename XML filename.
532  *
533  * Saves observations into a XML file. Please refer to the read() method for
534  * more information about the structure of the XML file.
535  ***************************************************************************/
536 void GObservations::save(const GFilename& filename) const
537 {
538  // Declare empty XML document
539  GXml xml;
540 
541  // Write observations into XML file
542  write(xml);
543 
544  // Save XML document
545  xml.save(filename);
546 
547  // Return
548  return;
549 }
550 
551 
552 /***********************************************************************//**
553  * @brief Read observations from XML document
554  *
555  * @param[in] xml XML document.
556  *
557  * @exception GException::invalid_value
558  * Invalid instrument encountered in XML file.
559  *
560  * Reads observations from the first observation list that is found in the
561  * XML document. The decoding of the instrument specific observation
562  * definition is done within the observation's GObservation::read() method.
563  * The following file structure is expected:
564  *
565  * <observation_list title="observation library">
566  * <observation name="..." id="..." instrument="...">
567  * ...
568  * </observation>
569  * <observation name="..." id="..." instrument="...">
570  * ...
571  * </observation>
572  * ...
573  * </observation_list>
574  *
575  * The @p name and @p id attributes allow for a unique identification of an
576  * observation within the observation container. The @p instrument
577  * attributes specifies the instrument to which the observation applies.
578  * This attribute will be used to allocate the appropriate instrument
579  * specific GObservation class variant by using the GObservationRegistry
580  * class.
581  *
582  * The structure within the @p observation tag is defined by the instrument
583  * specific GObservation class.
584  *
585  * @todo Observation names and IDs are not verified so far for uniqueness.
586  * This would be required to achieve an unambiguous update of parameters
587  * in an already existing XML file when using the write method.
588  ***************************************************************************/
589 void GObservations::read(const GXml& xml)
590 {
591  // Get pointer on observation library
592  const GXmlElement* lib = xml.element("observation_list", 0);
593 
594  // Loop over all observations
595  int n = lib->elements("observation");
596  for (int i = 0; i < n; ++i) {
597 
598  // Get pointer on observation
599  const GXmlElement* obs = lib->element("observation", i);
600 
601  // Get attributes
602  std::string name = obs->attribute("name");
603  std::string id = obs->attribute("id");
604  std::string instrument = obs->attribute("instrument");
605 
606  // Allocate instrument specific observation
607  GObservationRegistry registry;
608  GObservation* ptr = registry.alloc(instrument);
609 
610  // If observation is valid then read its definition from XML file
611  if (ptr != NULL) {
612 
613  // Read definition
614  ptr->read(*obs);
615 
616  // Set attributes
617  ptr->name(name);
618  ptr->id(id);
619 
620  } // endif: observation was valid
621 
622  // ... otherwise throw an exception
623  else {
624  std::string msg = "Instrument \""+instrument+"\" unknown. The "
625  "following instruments are available: "+
626  registry.content()+". Please specify one of "
627  "the available instruments.";
628  throw GException::invalid_value(G_READ, msg);
629  }
630 
631  // Append observation to container
632  append(*ptr);
633 
634  // Free observation (the append method made a deep copy)
635  delete ptr;
636 
637  } // endfor: looped over all observations
638 
639  // Return
640  return;
641 }
642 
643 
644 /***********************************************************************//**
645  * @brief Write observations into XML document
646  *
647  * @param[in] xml XML document.
648  *
649  * Write observations into the first observation library that is found in the
650  * XML document. In case that no observation library exists, one is added to
651  * the document. Please refer to the read() method for more information
652  * about the structure of the XML document.
653  ***************************************************************************/
654 void GObservations::write(GXml& xml) const
655 {
656  // If there is no observation library then append one
657  if (xml.elements("observation_list") == 0) {
658  xml.append(GXmlElement("observation_list title=\"observation list\""));
659  }
660 
661  // Get pointer on observation library
662  GXmlElement* lib = xml.element("observation_list", 0);
663 
664  // Write all observations into library
665  for (int i = 0; i < size(); ++i) {
666 
667  // Initialise pointer on observation
668  GXmlElement* obs = NULL;
669 
670  // Search corresponding observation
671  int n = xml.elements("observation");
672  for (int k = 0; k < n; ++k) {
673  GXmlElement* element = xml.element("observation", k);
674  if (element->attribute("name") == m_obs[i]->name() &&
675  element->attribute("id") == m_obs[i]->id() &&
676  element->attribute("instrument") == m_obs[i]->instrument()) {
677  obs = element;
678  break;
679  }
680  }
681 
682  // If no observation with corresponding name, ID and instrument was
683  // found then append one now
684  if (obs == NULL) {
685  obs = lib->append("observation");
686  obs->attribute("name", m_obs[i]->name());
687  obs->attribute("id", m_obs[i]->id());
688  obs->attribute("instrument", m_obs[i]->instrument());
689  }
690 
691  // Write now observation
692  m_obs[i]->write(*obs);
693 
694  } // endfor: looped over all observaitons
695 
696  // Return
697  return;
698 }
699 
700 
701 /***********************************************************************//**
702  * @brief Load models from XML file
703  *
704  * @param[in] filename XML filename.
705  *
706  * Loads the models from a XML file. Please refer to the GModels::read()
707  * method for more information about the expected structure of the XML
708  * file.
709  ***************************************************************************/
710 void GObservations::models(const GFilename& filename)
711 {
712  // Load models
713  m_models.load(filename);
714 
715  // Return
716  return;
717 }
718 
719 
720 /***********************************************************************//**
721  * @brief Optimize model parameters using optimizer
722  *
723  * @param[in] opt Optimizer.
724  *
725  * Optimizes the free parameters of the models by using the optimizer
726  * that has been provided by the @p opt argument.
727  ***************************************************************************/
729 {
730  // Extract optimizer parameter container from model container
731  GOptimizerPars pars = m_models.pars();
732 
733  // Optimize model parameters
734  opt.optimize(m_fct, pars);
735 
736  // Return
737  return;
738 }
739 
740 
741 /***********************************************************************//**
742  * @brief Computes parameter errors using optimizer
743  *
744  * @param[in] opt Optimizer.
745  *
746  * Calculates the errors of the free parameters of the models by using
747  * the optimizer that has been provided by the @p opt argument.
748  ***************************************************************************/
750 {
751  // Extract optimizer parameter container from model container
752  GOptimizerPars pars = m_models.pars();
753 
754  // Compute model parameter errors
755  opt.errors(m_fct, pars);
756 
757  // Return
758  return;
759 }
760 
761 
762 /***********************************************************************//**
763  * @brief Computes parameter errors using hessian matrix and optimizer
764  *
765  * Calculates the errors of the free model parameters as the square roots
766  * of the diagonal elements of the inverse Hessian matrix. The Hessian
767  * matrix is computed using finite differences.
768  ***************************************************************************/
770 {
771  // Extract optimizer parameter container from model container
772  GOptimizerPars pars = m_models.pars();
773 
774  // Get number of parameters
775  int npars = pars.size();
776 
777  // Compute Hessian matrix
778  GMatrixSparse hessian = m_fct.hessian(pars);
779 
780  // Signal no diagonal element loading
781  bool diag_loaded = false;
782 
783  // Loop over error computation (maximum 2 passes)
784  for (int i = 0; i < 2; ++i) {
785 
786  // Solve: hessian * X = unit
787  try {
788  GMatrixSparse decomposition = hessian.cholesky_decompose(true);
789  GVector unit(npars);
790  for (int ipar = 0; ipar < npars; ++ipar) {
791  unit[ipar] = 1.0;
792  GVector x = decomposition.cholesky_solver(unit, true);
793  if (x[ipar] >= 0.0) {
794  pars[ipar]->factor_error(sqrt(x[ipar]));
795  }
796  else {
797  pars[ipar]->factor_error(0.0);
798  }
799  unit[ipar] = 0.0;
800  }
801  }
802  catch (GException::runtime_error &e) {
803 
804  // Load diagonal if this has not yet been tried
805  if (!diag_loaded) {
806 
807  // Flag errors as inaccurate
808  std::cout << "Non-Positive definite hessian matrix encountered."
809  << std::endl;
810  std::cout << "Load diagonal elements with 1e-10."
811  << " Fit errors may be inaccurate."
812  << std::endl;
813 
814  // Try now with diagonal loaded matrix
815  for (int ipar = 0; ipar < npars; ++ipar) {
816  hessian(ipar,ipar) += 1.0e-10;
817  }
818 
819  // Signal loading
820  diag_loaded = true;
821 
822  // Try again
823  continue;
824 
825  } // endif: diagonal has not yet been loaded
826 
827  // ... otherwise signal an error
828  else {
829  std::cout << "Non-Positive definite hessian matrix encountered,"
830  << " even after diagonal loading." << std::endl;
831  break;
832  }
833  }
834  catch (std::exception &e) {
835  throw;
836  }
837 
838  // If no error occured then break now
839  break;
840 
841  } // endfor: looped over error computation
842 
843  // Return
844  return;
845 }
846 
847 
848 /***********************************************************************//**
849  * @brief Evaluate function
850  *
851  * Evaluates the likelihood funtion at the actual set of parameters.
852  ***************************************************************************/
854 {
855  // Extract optimizer parameter container from model container
856  GOptimizerPars pars = m_models.pars();
857 
858  // Compute log-likelihood
859  m_fct.eval(pars);
860 
861  // Return
862  return;
863 }
864 
865 
866 /***********************************************************************//**
867  * @brief Return total number of observed events
868  *
869  * @return Total number of observed events.
870  *
871  * Returns the total number of observed events that is container in the
872  * observation container.
873  ***************************************************************************/
875 {
876  // Initialise number of observed events
877  int nobserved = 0;
878 
879  // Compute number of observed events
880  for (int i = 0; i < size(); ++i) {
881  nobserved += (*this)[i]->nobserved();
882  }
883 
884  // Return number of observed events
885  return nobserved;
886 }
887 
888 
889 /***********************************************************************//**
890  * @brief Return number of predicted events for model with a given @p name
891  *
892  * @param[in] name Model name.
893  * @param[in] instrument Instrument name.
894  * @return Total number of predicted events for model.
895  *
896  * @exception GException::invalid_argument
897  * Model with @p name not found.
898  *
899  * Returns the total number of predicted events for the model with a given
900  * @p name.
901  ***************************************************************************/
902 double GObservations::npred(const std::string& name,
903  const std::string& instrument) const
904 {
905  // Initialise number of predicted events
906  double npred = 0.0;
907 
908  // Throw an exception if model does not exist
909  if (!m_models.contains(name)) {
910  std::string msg = "Model \""+name+"\" not found in observation "
911  "container. Please specify a valid model name.";
913  }
914 
915  // Get pointer to model
916  const GModel* model = m_models[name];
917 
918  // Continue only if pointer is valid
919  if (model != NULL) {
920  npred = this->npred(*model, instrument);
921  }
922 
923  // Return number of predicted events
924  return npred;
925 }
926 
927 
928 /***********************************************************************//**
929  * @brief Return number of predicted events for a given model
930  *
931  * @param[in] model Model.
932  * @param[in] instrument Instrument name.
933  * @return Total number of predicted events for model.
934  *
935  * @exception GException::invalid_argument
936  * Model with @p name not found.
937  *
938  * Returns the total number of predicted events for the model with a given
939  * @p name.
940  ***************************************************************************/
941 double GObservations::npred(const GModel& model,
942  const std::string& instrument) const
943 {
944  // Initialise number of predicted events
945  double npred = 0.0;
946 
947  // Compute number of predicted events
948  for (int i = 0; i < size(); ++i) {
949 
950  // Skip observation if it does not correspond to the specified
951  // instrument
952  if (!instrument.empty()) {
953  if ((*this)[i]->instrument() != instrument) {
954  continue;
955  }
956  }
957 
958  // Add only relevant events
959  if (model.is_valid(instrument, "")) {
960  npred += (*this)[i]->npred(model);
961  }
962 
963  } // endfor: looped over observations
964 
965  // Return number of predicted events
966  return npred;
967 }
968 
969 
970 /***********************************************************************//**
971  * @brief Remove response cache for all models
972  *
973  * Remove response cache for all models in the container.
974  ***************************************************************************/
976 {
977  // Loop over all models and remove response cache for each model
978  for (int i = 0; i < m_models.size(); ++i) {
979  remove_response_cache(m_models[i]->name());
980  }
981 
982  // Return
983  return;
984 }
985 
986 
987 /***********************************************************************//**
988  * @brief Remove response cache for model
989  *
990  * @param[in] name Model name.
991  *
992  * Remove response cache for model @p name in the container.
993  ***************************************************************************/
994 void GObservations::remove_response_cache(const std::string& name)
995 {
996  // Loop over all observations and remove response cache for the
997  // specified source
998  for (int i = 0; i < size(); ++i) {
999  (*this)[i]->remove_response_cache(name);
1000  }
1001 
1002  // Return
1003  return;
1004 }
1005 
1006 
1007 /***********************************************************************//**
1008  * @brief Print observation list information
1009  *
1010  * @param[in] chatter Chattiness.
1011  * @return String containing observation list information
1012  ***************************************************************************/
1013 std::string GObservations::print(const GChatter& chatter) const
1014 {
1015  // Initialise result string
1016  std::string result;
1017 
1018  // Continue only if chatter is not silent
1019  if (chatter != SILENT) {
1020 
1021  // Append header
1022  result.append("=== GObservations ===");
1023 
1024  // Append information
1025  result.append("\n"+gammalib::parformat("Number of observations"));
1026  result.append(gammalib::str(size()));
1027  result.append("\n"+gammalib::parformat("Number of models"));
1028  result.append(gammalib::str(m_models.size()));
1029  result.append("\n"+gammalib::parformat("Number of observed events"));
1030  result.append(gammalib::str(nobserved()));
1031  result.append("\n"+gammalib::parformat("Number of predicted events"));
1032  result.append(gammalib::str(npred()));
1033 
1034  // Get reduced chatter level for detailed information
1035  GChatter reduced_chatter = gammalib::reduce(chatter);
1036 
1037  // Append observations for chatter >= EXPLICIT
1038  if (chatter >= EXPLICIT) {
1039  for (int i = 0; i < size(); ++i) {
1040  result.append("\n");
1041  result.append((*this)[i]->print(reduced_chatter));
1042  }
1043  }
1044 
1045  // Append models for chatter >= VERBOSE
1046  if (chatter >= VERBOSE) {
1047  result.append("\n"+m_models.print(reduced_chatter));
1048  }
1049 
1050  } // endif: chatter was not silent
1051 
1052  // Return result
1053  return result;
1054 }
1055 
1056 
1057 /*==========================================================================
1058  = =
1059  = Private methods =
1060  = =
1061  ==========================================================================*/
1062 
1063 /***********************************************************************//**
1064  * @brief Initialise class members
1065  ***************************************************************************/
1067 {
1068  // Initialise members
1069  m_obs.clear();
1070  m_models.clear();
1071  m_fct.set(this); //!< Makes sure that optimizer points to this instance
1072 
1073  // Return
1074  return;
1075 }
1076 
1077 
1078 /***********************************************************************//**
1079  * @brief Copy class members
1080  *
1081  * @param[in] obs Observation container.
1082  *
1083  * Copies all class members. Deep copies are created for the observations,
1084  * which is what is needed to have a real physical copy of the members.
1085  ***************************************************************************/
1087 {
1088  // Copy attributes
1089  m_models = obs.m_models;
1090  m_fct = obs.m_fct;
1091 
1092  // Copy observations
1093  m_obs.clear();
1094  for (int i = 0; i < obs.m_obs.size(); ++i) {
1095  m_obs.push_back((obs.m_obs[i]->clone()));
1096  }
1097 
1098  // Set likelihood function pointer to this instance. This makes sure
1099  // that the optimizer points to this instance
1100  m_fct.set(this);
1101 
1102  // Return
1103  return;
1104 }
1105 
1106 
1107 /***********************************************************************//**
1108  * @brief Delete class members
1109  *
1110  * Deallocates all observations. Since container classes that hold pointers
1111  * need to handle the proper deallocation of memory, we loop here over all
1112  * pointers and make sure that we deallocate all observations in the
1113  * container.
1114  ***************************************************************************/
1116 {
1117  // Free observations
1118  for (int i = 0; i < m_obs.size(); ++i) {
1119  delete m_obs[i];
1120  m_obs[i] = NULL;
1121  }
1122 
1123  // Return
1124  return;
1125 }
1126 
1127 
1128 /***********************************************************************//**
1129  * @brief Return observation index by instrument and identifier
1130  *
1131  * @param[in] instrument Instrument.
1132  * @param[in] id Observation identifier.
1133  * @return Observation index (-1 of instrument and observation identifier
1134  * has not been found)
1135  *
1136  * Returns observation index based on the specified @p instrument and
1137  * observation identifier @p id. If no observation with the specified
1138  * attributes has been found, the method returns -1.
1139  ***************************************************************************/
1140 int GObservations::get_index(const std::string& instrument,
1141  const std::string& id) const
1142 {
1143  // Initialise index
1144  int index = -1;
1145 
1146  // Search observation with specified instrument and id
1147  for (int i = 0; i < size(); ++i) {
1148  if ((m_obs[i]->instrument() == instrument) &&
1149  (m_obs[i]->id() == id)) {
1150  index = i;
1151  break;
1152  }
1153  }
1154 
1155  // Return index
1156  return index;
1157 }
Abstract model class.
Definition: GModel.hpp:100
GObservations * clone(void) const
Clone observations.
#define G_EXTEND
Abstract model base class interface definition.
virtual void read(const GXmlElement &xml)=0
GObservation * alloc(const std::string &name) const
Allocate observation of given name.
Sparse matrix class interface definition.
void save(const GFilename &filename) const
Save XML document into file.
Definition: GXml.cpp:581
virtual std::string instrument(void) const =0
XML class interface definition.
#define G_REMOVE
void save(const GFilename &filename) const
Save observations into XML file.
#define G_NPRED
Optimizer parameter container class.
XML element node class.
Definition: GXmlElement.hpp:48
#define G_READ
GOptimizerPars pars(void) const
Return optimizer parameter container.
Definition: GModels.cpp:878
GMatrixSparse hessian(const GOptimizerPars &pars)
Compute Hessian matrix.
GObservation * append(const GObservation &obs)
Append observation to container.
virtual int elements(void) const
Return number of GXMLElement children of node.
Definition: GXmlNode.cpp:586
bool contains(const std::string &instrument, const std::string &id) const
Signals if observation exists.
Gammalib tools definition.
void errors(GOptimizer &opt)
Computes parameter errors using optimizer.
void id(const std::string &id)
Set observation identifier.
Abstract optimizer abstract base class interface definition.
#define G_AT
const GModels & models(void) const
Return model container.
bool is_valid(const std::string &instruments, const std::string &ids) const
Verifies if model is valid for a given instrument and identifier.
Definition: GModel.cpp:570
virtual ~GObservations(void)
Destructor.
std::string print(const GChatter &chatter=NORMAL) const
Print observation list information.
GObservation * set(const int &index, const GObservation &obs)
Set observation in container.
void clear(void)
Clear object.
Definition: GModels.cpp:233
std::string content(void) const
Return list of names in registry.
Definition: GRegistry.cpp:54
void copy_members(const GObservations &obs)
Copy class members.
std::vector< GObservation * > m_obs
List of observations.
void errors_hessian(void)
Computes parameter errors using hessian matrix and optimizer.
GVector sqrt(const GVector &vector)
Computes square root of vector elements.
Definition: GVector.cpp:1358
void set(GObservations *obs)
Set observation container.
GObservations(void)
Void constructor.
GXmlElement * element(const int &index)
Return pointer to child element.
Definition: GXml.cpp:419
bool is_empty(void) const
Signals if there are no observations in container.
int elements(void) const
Return number of child elements in XML document root.
Definition: GXml.cpp:384
GModels m_models
List of models.
virtual void optimize(GOptimizerFunction &fct, GOptimizerPars &pars)=0
const GXmlAttribute * attribute(const int &index) const
Return attribute.
XML class.
Definition: GXml.hpp:172
Filename class.
Definition: GFilename.hpp:62
GObservations & operator=(const GObservations &obs)
Assignment operator.
double npred(void) const
Return total number of predicted events in models.
Interface definition for the observation registry class.
int get_index(const std::string &instrument, const std::string &id) const
Return observation index by instrument and identifier.
void read(const GXml &xml)
Read observations from XML document.
Observation registry class definition.
void init_members(void)
Initialise class members.
Abstract optimizer abstract base class.
Definition: GOptimizer.hpp:54
virtual GObservation * clone(void) const =0
Clones object.
GChatter
Definition: GTypemaps.hpp:33
void optimize(GOptimizer &opt)
Optimize model parameters using optimizer.
Abstract observation base class.
int size(void) const
Return number of observations in container.
void eval(void)
Evaluate function.
Observation container class.
GObservation * insert(const int &index, const GObservation &obs)
Insert observation into container.
#define G_APPEND
virtual void eval(const GOptimizerPars &pars)
Evaluate log-likelihood function.
virtual void errors(GOptimizerFunction &fct, GOptimizerPars &pars)=0
void name(const std::string &name)
Set observation name.
std::string url(void) const
Return Uniform Resource Locator (URL)
Definition: GFilename.hpp:189
int nobserved(void) const
Return total number of observed events.
#define G_INSERT
virtual GXmlElement * element(const int &index)
Return pointer to GXMLElement child.
Definition: GXmlNode.cpp:640
void load(const GFilename &filename)
Load models from XML file.
Definition: GModels.cpp:688
void extend(const GObservations &obs)
Append observations from observation container.
void write(GXml &xml) const
Write observations into XML document.
bool contains(const std::string &name) const
Signals if model name exists.
Definition: GModels.cpp:670
GChatter reduce(const GChatter &chatter)
Reduce chattiness by one level.
Definition: GTypemaps.hpp:65
Sparse matrix class definition.
std::string print(const GChatter &chatter=NORMAL) const
Print models.
Definition: GModels.cpp:972
GVector cholesky_solver(const GVector &vector, const bool &compress=true) const
Cholesky solver.
Exception handler interface definition.
Vector class.
Definition: GVector.hpp:46
int size(void) const
Return number of parameters in container.
void load(const GFilename &filename)
Load observations from XML file.
void clear(void)
Clear observations.
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
void free_members(void)
Delete class members.
#define G_SET
void remove_response_cache(void)
Remove response cache for all models.
GObservations::likelihood m_fct
Optimizer function.
Observations container class interface definition.
GMatrixSparse cholesky_decompose(const bool &compress=true) const
Return Cholesky decomposition.
GObservation * at(const int &index)
Return pointer to observation.
int size(void) const
Return number of models in container.
Definition: GModels.hpp:259
void remove(const int &index)
Remove observation from container.
Filename class interface definition.
GXmlNode * append(const GXmlNode &node)
Append child node to XML document root.
Definition: GXml.cpp:279
std::string str(const unsigned short int &value)
Convert unsigned short integer value into string.
Definition: GTools.cpp:489
void reserve(const int &num)
Reserves space for observations in container.