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