GammaLib  1.7.0.dev
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Macros Pages
GObservation.cpp
Go to the documentation of this file.
1 /***************************************************************************
2  * GObservation.cpp - Abstract observation base class *
3  * ----------------------------------------------------------------------- *
4  * copyright (C) 2008-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 GObservation.cpp
23  * @brief Abstract observation base 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 "GIntegral.hpp"
34 #include "GDerivative.hpp"
35 #include "GVector.hpp"
36 #include "GMatrixSparse.hpp"
37 #include "GObservation.hpp"
38 #include "GResponse.hpp"
39 #include "GEventCube.hpp"
40 #include "GEventList.hpp"
41 #include "GEventBin.hpp"
42 #include "GModel.hpp"
43 #include "GModels.hpp"
44 #include "GModelPar.hpp"
45 #include "GModelSky.hpp"
46 #include "GModelData.hpp"
47 
48 /* __ Method name definitions ____________________________________________ */
49 #define G_LIKELIHOOD "GObservation::likelihood(GModels&, GVector*,"\
50  " GMatrixSparse*, double*)"
51 #define G_MODEL "GObservation::model(GModels&, GPointing&,"\
52  " GInstDir&, GEnergy&, GTime&, GVector*)"
53 #define G_EVENTS "GObservation::events()"
54 #define G_NPRED "GObservation::npred(GModel&)"
55 #define G_NPRED_SPEC "GObservation::npred_spec(GModel&, GTime&)"
56 
57 /* __ Constants __________________________________________________________ */
58 const double minmod = 1.0e-100; //!< Minimum model value
59 const double minerr = 1.0e-100; //!< Minimum statistical error
60 
61 /* __ Macros _____________________________________________________________ */
62 
63 /* __ Coding definitions _________________________________________________ */
64 #define G_LN_ENERGY_INT //!< ln(E) variable substitution for integration
65 
66 /* __ Debug definitions __________________________________________________ */
67 //#define G_OPT_DEBUG
68 
69 
70 /*==========================================================================
71  = =
72  = Constructors/destructors =
73  = =
74  ==========================================================================*/
75 
76 /***********************************************************************//**
77  * @brief Void constructor
78  ***************************************************************************/
80 {
81  // Initialise members
82  init_members();
83 
84  // Return
85  return;
86 }
87 
88 
89 /***********************************************************************//**
90  * @brief Copy constructor
91  *
92  * @param[in] obs Observation.
93  *
94  * Instantiate the class by copying from an existing observation.
95  ***************************************************************************/
97 {
98  // Initialise members
99  init_members();
100 
101  // Copy members
102  copy_members(obs);
103 
104  // Return
105  return;
106 }
107 
108 
109 /***********************************************************************//**
110  * @brief Destructor
111  ***************************************************************************/
113 {
114  // Free members
115  free_members();
116 
117  // Return
118  return;
119 }
120 
121 
122 /*==========================================================================
123  = =
124  = Operators =
125  = =
126  ==========================================================================*/
127 
128 /***********************************************************************//**
129  * @brief Assignment operator
130  *
131  * @param[in] obs Observation.
132  * @return Observation.
133  *
134  * Assign observation.
135  ***************************************************************************/
137 {
138  // Execute only if object is not identical
139  if (this != &obs) {
140 
141  // Free members
142  free_members();
143 
144  // Initialise members
145  init_members();
146 
147  // Copy members
148  copy_members(obs);
149 
150  } // endif: object was not identical
151 
152  // Return this object
153  return *this;
154 }
155 
156 
157 /*==========================================================================
158  = =
159  = Public methods =
160  = =
161  ==========================================================================*/
162 
163 /***********************************************************************//**
164  * @brief Compute likelihood function
165  *
166  * @param[in] models Models.
167  * @param[in,out] gradient Pointer to gradients.
168  * @param[in,out] curvature Pointer to curvature matrix.
169  * @param[in,out] npred Pointer to Npred value.
170  * @return Likelihood.
171  *
172  * Computes the likelihood for a specified set of models. The method also
173  * returns the gradients, the curvature matrix, and the number of events
174  * that are predicted by all models.
175  ***************************************************************************/
176 double GObservation::likelihood(const GModels& models,
177  GVector* gradient,
178  GMatrixSparse* curvature,
179  double* npred) const
180 {
181  // Initialise likelihood value
182  double value = 0.0;
183 
184  // Extract statistic for this observation
185  std::string statistic = gammalib::toupper(this->statistic());
186 
187  // Unbinned analysis
188  if (dynamic_cast<const GEventList*>(events()) != NULL) {
189 
190  // Poisson statistic
191  if ((statistic == "POISSON") || (statistic == "CSTAT")) {
192 
193  // Update the log-likelihood
194  value = likelihood_poisson_unbinned(models,
195  gradient,
196  curvature,
197  npred);
198 
199  } // endif: Poisson statistic
200 
201  // ... otherwise throw an exception
202  else {
203  std::string msg = "Invalid statistic \""+statistic+"\". Unbinned "
204  "optimization requires \"POISSON\" or \"CSTAT\" "
205  "statistic.";
207  }
208 
209  } // endif: unbinned analysis
210 
211  // ... or binned analysis
212  else {
213 
214  // Poisson statistic
215  if ((statistic == "POISSON") || (statistic == "CSTAT")) {
216  value = likelihood_poisson_binned(models,
217  gradient,
218  curvature,
219  npred);
220  }
221 
222  // ... or Gaussian statistic
223  else if ((statistic == "GAUSSIAN") || (statistic == "CHI2")) {
224  value = likelihood_gaussian_binned(models,
225  gradient,
226  curvature,
227  npred);
228  }
229 
230  // ... or unsupported
231  else {
232  std::string msg = "Invalid statistic \""+statistic+"\". Binned "
233  "optimization requires \"POISSON\", \"CSTAT\", "
234  "\"GAUSSIAN\" or \"CHI2\" statistic.";
236  }
237 
238  } // endelse: binned analysis
239 
240  // Return likelihood
241  return value;
242 }
243 
244 
245 /***********************************************************************//**
246  * @brief Return model value and (optionally) gradient
247  *
248  * @param[in] models Model descriptor.
249  * @param[in] event Observed event.
250  * @param[out] gradient Pointer to gradient vector (optional).
251  *
252  * @exception GException::invalid_value
253  * Dimension of gradient vector mismatches number of parameters.
254  *
255  * Implements generic model and gradient evaluation for an observation. The
256  * model gives the probability for an event to occur with a given instrument
257  * direction, at a given energy and at a given time. The gradient is the
258  * parameter derivative of this probability. If NULL is passed for the
259  * gradient vector, then gradients will not be computed.
260  *
261  * The method will only operate on models for which the list of instruments
262  * and observation identifiers matches those of the observation. Models that
263  * do not match will be skipped.
264  ***************************************************************************/
265 double GObservation::model(const GModels& models,
266  const GEvent& event,
267  GVector* gradient) const
268 {
269  // Initialise method variables
270  double model = 0.0; // Reset model value
271  int igrad = 0; // Reset gradient counter
272  int grad_size = 0; // Reset gradient size
273 
274  // If gradient is available then reset gradient vector elements to 0
275  // and determine vector size
276  if (gradient != NULL) {
277  (*gradient) = 0.0;
278  grad_size = gradient->size();
279  }
280 
281  // Loop over models
282  for (int i = 0; i < models.size(); ++i) {
283 
284  // Get model pointer. Continue only if pointer is valid
285  const GModel* mptr = models[i];
286  if (mptr != NULL) {
287 
288  // Continue only if model applies to specific instrument and
289  // observation identifier
290  if (mptr->is_valid(instrument(), id())) {
291 
292  // Compute value and add to model. If energy dispersion
293  // is used, don't compute model gradients as we cannot
294  // use them. This is somehow a kluge, but makes the
295  // code faster
296  model += mptr->eval(event, *this, !response()->use_edisp());
297 
298  // Optionally determine model gradients. If the model has a
299  // gradient then use it, unless we have energy dispersion.
300  // For energy dispersion, no gradients are available as we
301  // have not implemented code that integrates the gradients
302  // over the energy dispersion. Here it's simpler to just use
303  // numerical gradients.
304  if (gradient != NULL) {
305  for (int ipar = 0; ipar < mptr->size(); ++ipar) {
306 
307  // Get reference to model parameter
308  const GModelPar& par = (*mptr)[ipar];
309 
310  // Make sure that we have a slot for the gradient
311  #if defined(G_RANGE_CHECK)
312  if (igrad >= grad_size) {
313  std::string msg = "Vector has not enough elements "
314  "to store the model parameter "
315  "gradients. "+
316  gammalib::str(models.npars())+
317  " elements requested while vector "
318  "only contains "+
319  gammalib::str(gradient->size())+
320  " elements.";
322  }
323  #endif
324 
325  if (par.is_free()) {
326  if (par.has_grad() && !response()->use_edisp()) {
327  (*gradient)[igrad+ipar] = par.factor_gradient();
328  }
329  else {
330  (*gradient)[igrad+ipar] = model_grad(*mptr, par, event);
331  }
332  }
333  else {
334  (*gradient)[igrad+ipar] = 0.0;
335  }
336  }
337  }
338 
339  } // endif: model component was valid for instrument
340 
341  // Increment parameter counter for gradients
342  igrad += mptr->size();
343 
344  } // endif: model was valid
345 
346  } // endfor: Looped over models
347 
348  // Return
349  return model;
350 }
351 
352 
353 /***********************************************************************//**
354  * @brief Return total number of observed events
355  *
356  * @returns Total number of observed events.
357  *
358  * Returns the total number of observed events. If the observation does not
359  * contain any events the method returns zero.
360  ***************************************************************************/
361 int GObservation::nobserved(void) const
362 {
363  // Initialise number of observed events
364  int nobserved = 0;
365 
366  // Extract number of observed events
367  if (m_events != NULL) {
368  nobserved = m_events->number();
369  }
370 
371  // Return number of observed events
372  return nobserved;
373 }
374 
375 
376 /***********************************************************************//**
377  * @brief Return total number (and optionally gradient) of predicted counts
378  * for all models
379  *
380  * @param[in] models Models.
381  * @param[out] gradient Model parameter gradients (optional).
382  *
383  * @exception GException::gradient_par_mismatch
384  * Dimension of gradient vector mismatches number of parameters.
385  *
386  * Returns the total number of predicted counts within the analysis region.
387  * If NULL is passed for the gradient vector then gradients will not be
388  * computed.
389  *
390  * The method will only operate on models for which the list of instruments
391  * and observation identifiers matches those of the observation. Models that
392  * do not match will be skipped.
393  ***************************************************************************/
394 double GObservation::npred(const GModels& models, GVector* gradient) const
395 {
396  // Verify that gradient vector and models have the same dimension
397  #if defined(G_RANGE_CHECK)
398  if (gradient != NULL) {
399  if (models.npars() != gradient->size()) {
401  gradient->size(),
402  models.npars());
403  }
404  }
405  #endif
406 
407  // Initialise
408  double npred = 0.0; // Reset predicted number of counts
409  int igrad = 0; // Reset gradient counter
410 
411  // If gradient is available then reset gradient vector elements to 0
412  if (gradient != NULL) {
413  (*gradient) = 0.0;
414  }
415 
416  // Loop over models
417  for (int i = 0; i < models.size(); ++i) {
418 
419  // Get model pointer. Continue only if pointer is valid
420  const GModel* mptr = models[i];
421  if (mptr != NULL) {
422 
423  // Continue only if model applies to specific instrument and
424  // observation identifier
425  if (mptr->is_valid(instrument(), id())) {
426 
427  // Determine Npred for model
428  npred += this->npred(*mptr);
429 
430  // Optionally determine Npred gradients
431  if (gradient != NULL) {
432  for (int k = 0; k < mptr->size(); ++k) {
433  const GModelPar& par = (*mptr)[k];
434  (*gradient)[igrad+k] = npred_grad(*mptr, par);
435  }
436  }
437 
438  } // endif: model component was valid for instrument
439 
440  // Increment parameter counter for gradient
441  igrad += mptr->size();
442 
443  } // endif: model was valid
444 
445  } // endfor: Looped over models
446 
447  // Return prediction
448  return npred;
449 }
450 
451 
452 /***********************************************************************//**
453  * @brief Return total number of predicted counts for one model
454  *
455  * @param[in] model Gamma-ray source model.
456  *
457  * @exception GException::gti_invalid
458  * Good Time Interval is invalid.
459  *
460  * Computes
461  *
462  * \f[
463  * N_{\rm pred} = \int_{\rm GTI} \int_{E_{\rm bounds}} \int_{\rm ROI}
464  * P(p',E',t') \, dp' \, dE' \, dt'
465  * \f]
466  *
467  * of the event probability
468  *
469  * \f[
470  * P(p',E',t') = \int \int \int
471  * S(p,E,t) \times R(p',E',t'|p,E,t) \, dp \, dE \, dt
472  * \f]
473  *
474  * where
475  * \f$S(p,E,t)\f$ is the source model,
476  * \f$R(p',E',t'|p,E,t)\f$ is the instrument response function,
477  * \f$p'\f$ is the measured photon direction,
478  * \f$E'\f$ is the measured photon energy,
479  * \f$t'\f$ is the measured photon arrival time,
480  * \f$p\f$ is the true photon arrival direction,
481  * \f$E\f$ is the true photon energy, and
482  * \f$t\f$ is the true photon arrival time.
483  *
484  * This method performs the time integration over the Good Time Intervals
485  * \f${\rm GTI}\f$. Note that MET is used for the time integration interval.
486  * This, however, is no specialisation since npred_grad_kern_spec::eval()
487  * converts the argument back in a GTime object by assuming that the argument
488  * is in MET, hence the correct time system will be used at the end by the
489  * method.
490  ***************************************************************************/
491 double GObservation::npred(const GModel& model) const
492 {
493  // Initialise result
494  double npred = 0.0;
495 
496  // Continue only if model applies to specific instrument and
497  // observation identifier
498  if (model.is_valid(instrument(), id())) {
499 
500  // Case A: If the model is constant then integrate analytically
501  if (model.is_constant()) {
502 
503  // Evaluate model at first start time and multiply by ontime
504  double ontime = events()->gti().ontime();
505 
506  // Integrate only if ontime is positive
507  if (ontime > 0.0) {
508 
509  // Integration is a simple multiplication by the time
510  npred = npred_spec(model, events()->gti().tstart()) * ontime;
511 
512  }
513 
514  } // endif: model was constant
515 
516  // ... otherwise integrate temporally
517  else {
518 
519  // Loop over GTIs
520  for (int i = 0; i < events()->gti().size(); ++i) {
521 
522  // Set integration interval in seconds
523  double tstart = events()->gti().tstart(i).secs();
524  double tstop = events()->gti().tstop(i).secs();
525 
526  // Throw exception if time interval is not valid
527  if (tstop <= tstart) {
529  events()->gti().tstart(i),
530  events()->gti().tstop(i));
531  }
532 
533  // Setup integration function
534  GObservation::npred_kern integrand(this, &model);
535  GIntegral integral(&integrand);
536 
537  // Do Romberg integration
538  npred += integral.romberg(tstart, tstop);
539 
540  } // endfor: looped over GTIs
541 
542  } // endelse: integrated temporally
543 
544  } // endif: model was valid
545 
546  // Return Npred
547  return npred;
548 }
549 
550 
551 /***********************************************************************//**
552  * @brief Returns parameter gradient of model for a given event
553  *
554  * @param[in] model Model.
555  * @param[in] par Model parameter.
556  * @param[in] event Event.
557  *
558  * This method uses a robust but simple difference method to estimate
559  * parameter gradients that have not been provided by the model. We use here
560  * a simple method as this method is likely used for spatial model parameters,
561  * and the spatial model may eventually be noisy due to numerical integration
562  * limits.
563  *
564  * The step size for the simple method has been fixed to 0.0002, which
565  * corresponds to about 1 arcsec for parameters that are given in degrees.
566  * The reasoning behind this value is that parameters that use numerical
567  * gradients are typically angles, such as for example the position, and
568  * we want to achieve arcsec precision with this method.
569  ***************************************************************************/
570 double GObservation::model_grad(const GModel& model,
571  const GModelPar& par,
572  const GEvent& event) const
573 {
574  // Initialise gradient
575  double grad = 0.0;
576 
577  // Compute gradient only if parameter is free
578  if (par.is_free()) {
579 
580  // Get non-const model pointer
581  GModelPar* ptr = const_cast<GModelPar*>(&par);
582 
583  // Save current model parameter
584  GModelPar current = par;
585 
586  // Get actual parameter value
587  double x = par.factor_value();
588 
589  // Set fixed step size for computation of derivative.
590  // By default, the step size is fixed to 0.0002.
591  const double step_size = 0.0002; // ~1 arcsec
592  double h = step_size;
593 
594  // Re-adjust the step-size h in case that the initial step size is
595  // larger than the allowed parameter range
596  if (par.has_min() && par.has_max()) {
597  double par_h = par.factor_max() - par.factor_min();
598  if (par_h < h) {
599  h = par_h;
600  }
601  }
602 
603  // Continue only if step size is positive
604  if (h > 0.0) {
605 
606  // Remove any boundaries to avoid limitations
607  //ptr->remove_range(); // Not needed in principle
608 
609  // Setup derivative function
610  GObservation::model_func function(this, &model, ptr, &event);
611  GDerivative derivative(&function);
612 
613  // If we are too close to the minimum boundary use a right sided
614  // difference ...
615  if (par.has_min() && ((x-par.factor_min()) < h)) {
616  grad = derivative.right_difference(x, h);
617  }
618 
619  // ... otherwise if we are too close to the maximum boundary use
620  // a left sided difference ...
621  else if (par.has_max() && ((par.factor_max()-x) < h)) {
622  grad = derivative.left_difference(x, h);
623  }
624 
625  // ... otherwise use a symmetric difference
626  else {
627  grad = derivative.difference(x, h);
628  }
629 
630  } // endif: step size was positive
631 
632  // Restore current model parameter
633  *ptr = current;
634 
635  } // endif: model parameter was free
636 
637  // Return gradient
638  return grad;
639 }
640 
641 
642 /***********************************************************************//**
643  * @brief Returns parameter gradient of Npred
644  *
645  * @param[in] model Gamma-ray source model.
646  * @param[in] par Model parameter.
647  *
648  * Computes
649  *
650  * \f[
651  * \frac{{\rm d} N_{\rm pred}}{{\rm d} a_i}
652  * \f]
653  *
654  * where
655  *
656  * \f[
657  * N_{\rm pred} = \int_{\rm GTI} \int_{E_{\rm bounds}} \int_{\rm ROI}
658  * P(p',E',t') \, dp' \, dE' \, dt'
659  * \f]
660  *
661  * is the integral of the event probability
662  *
663  * \f[
664  * P(p',E',t') = \int \int \int
665  * S(p,E,t) \times R(p',E',t'|p,E,t) \, dp \, dE \, dt
666  * \f]
667  *
668  * and
669  * \f$a_i\f$ is the model parameter \f$i\f$.
670  * Furthermore
671  * \f$S(p,E,t)\f$ is the source model,
672  * \f$R(p',E',t'|p,E,t)\f$ is the instrument response function,
673  * \f$p'\f$ is the measured photon direction,
674  * \f$E'\f$ is the measured photon energy,
675  * \f$t'\f$ is the measured photon arrival time,
676  * \f$p\f$ is the true photon arrival direction,
677  * \f$E\f$ is the true photon energy, and
678  * \f$t\f$ is the true photon arrival time.
679  *
680  * This method uses a robust but simple difference method to estimate
681  * parameter gradients that have not been provided by the model. We use here
682  * a simple method as this method is likely used for spatial model parameters,
683  * and the spatial model may eventually be noisy due to numerical integration
684  * limits.
685  *
686  * The step size for the simple method has been fixed to 0.0002.
687  ***************************************************************************/
688 double GObservation::npred_grad(const GModel& model, const GModelPar& par) const
689 {
690  // Initialise result
691  double grad = 0.0;
692 
693  // Compute gradient only if parameter is free
694  if (par.is_free()) {
695 
696  // Get non-const model pointer
697  GModelPar* ptr = const_cast<GModelPar*>(&par);
698 
699  // Save current model parameter
700  GModelPar current = par;
701 
702  // Get actual parameter value
703  double x = par.factor_value();
704 
705  // Set fixed step size to 0.0002 for computation of derivative.
706  const double step_size = 0.0002;
707  double h = step_size;
708 
709  // Re-adjust the step-size h in case that the initial step size is
710  // larger than the allowed parameter range
711  if (par.has_min() && par.has_max()) {
712  double par_h = par.factor_max() - par.factor_min();
713  if (par_h < h) {
714  h = par_h;
715  }
716  }
717 
718  // Continue only if step size is positive
719  if (h > 0.0) {
720 
721  // Remove any boundaries to avoid limitations
722  //ptr->remove_range(); // Not needed in principle
723 
724  // Setup derivative function
725  GObservation::npred_func function(this, &model, ptr);
726  GDerivative derivative(&function);
727 
728  // If we are too close to the minimum boundary use a right sided
729  // difference ...
730  if (par.has_min() && ((x-par.factor_min()) < h)) {
731  grad = derivative.right_difference(x, h);
732  }
733 
734  // ... otherwise if we are too close to the maximum boundary use
735  // a left sided difference ...
736  else if (par.has_max() && ((par.factor_max()-x) < h)) {
737  grad = derivative.left_difference(x, h);
738  }
739 
740  // ... otherwise use a symmetric difference
741  else {
742  grad = derivative.difference(x, h);
743  }
744 
745  } // endif: step size was positive
746 
747  // Restore current model parameter
748  *ptr = current;
749 
750  } // endif: model parameter was free
751 
752  // Return result
753  return grad;
754 }
755 
756 
757 /***********************************************************************//**
758  * @brief Response cache removal hook
759  *
760  * @param[in] name Model name.
761  *
762  * This is a response cache removal hook that is called by the
763  * GObservations::remove_response_cache() methods and that can be used to
764  * remove the model @p name from a response cache, if implemented in the
765  * derived class.
766  ***************************************************************************/
767 void GObservation::remove_response_cache(const std::string& name)
768 {
769  // Return
770  return;
771 }
772 
773 
774 /***********************************************************************//**
775  * @brief Set event container
776  *
777  * @param[in] events Event container.
778  *
779  * Set the event container for this observation by cloning the container
780  * specified in the argument.
781  ***************************************************************************/
782 void GObservation::events(const GEvents& events)
783 {
784  // Free existing event container only if it differs from current event
785  // container. This prevents unintential deallocation of the argument
786  if ((m_events != NULL) && (m_events != &events)) {
787  delete m_events;
788  }
789 
790  // Clone events
791  m_events = events.clone();
792 
793  // Return
794  return;
795 }
796 
797 
798 /***********************************************************************//**
799  * @brief Return events
800  *
801  * @exception GException::no_events
802  * No events allocated for observation.
803  *
804  * Returns pointer to events.
805  ***************************************************************************/
807 {
808  // Throw an exception if the event container is not valid
809  if (m_events == NULL) {
810  std::string msg = "No events allocated for observation.";
812  }
813 
814  // Return pointer to events
815  return (m_events);
816 }
817 
818 
819 /***********************************************************************//**
820  * @brief Return events (const version)
821  *
822  * @exception GException::no_events
823  * No events allocated for observation.
824  *
825  * Returns const pointer to events.
826  ***************************************************************************/
827 const GEvents* GObservation::events(void) const
828 {
829  // Throw an exception if the event container is not valid
830  if (m_events == NULL) {
831  std::string msg = "No events allocated for observation.";
833  }
834 
835  // Return pointer to events
836  return (m_events);
837 }
838 
839 
840 /*==========================================================================
841  = =
842  = Private methods =
843  = =
844  ==========================================================================*/
845 
846 /***********************************************************************//**
847  * @brief Initialise class members
848  ***************************************************************************/
850 {
851  // Initialise members
852  m_name.clear();
853  m_id.clear();
854  m_statistic = "cstat";
855  m_events = NULL;
856 
857  // Return
858  return;
859 }
860 
861 
862 /***********************************************************************//**
863  * @brief Copy class members
864  *
865  * @param[in] obs Observation.
866  *
867  * Copy members from an observation.
868  ***************************************************************************/
870 {
871  // Copy members
872  m_name = obs.m_name;
873  m_id = obs.m_id;
874  m_statistic = obs.m_statistic;
875 
876  // Clone members
877  m_events = (obs.m_events != NULL) ? obs.m_events->clone() : NULL;
878 
879  // Return
880  return;
881 }
882 
883 
884 /***********************************************************************//**
885  * @brief Delete class members
886  ***************************************************************************/
888 {
889  // Free members
890  if (m_events != NULL) delete m_events;
891 
892  // Signal free pointers
893  m_events = NULL;
894 
895  // Return
896  return;
897 }
898 
899 
900 /*==========================================================================
901  = =
902  = Likelihood methods =
903  = =
904  ==========================================================================*/
905 
906 /***********************************************************************//**
907  * @brief Evaluate log-likelihood function for Poisson statistic and
908  * unbinned analysis (version with working arrays)
909  *
910  * @param[in] models Models.
911  * @param[in,out] gradient Gradient.
912  * @param[in,out] curvature Curvature matrix.
913  * @param[in,out] npred Number of predicted events.
914  * @return Likelihood value.
915  *
916  * This method evaluates the -(log-likelihood) function for parameter
917  * optimisation using unbinned analysis and Poisson statistic.
918  * The -(log-likelihood) function is given by
919  *
920  * \f[
921  * L = N_{\rm pred} - \sum_i \log e_i
922  * \f]
923  *
924  * where
925  * \f$N_{\rm pred}\f$ is the number of events predicted by the model, and
926  * the sum is taken over all events. This method also computes the
927  * parameter gradients
928  * \f$\delta L/dp\f$
929  * and the curvature matrix
930  * \f$\delta^2 L/dp_1 dp_2\f$.
931  ***************************************************************************/
933  GVector* gradient,
934  GMatrixSparse* curvature,
935  double* npred) const
936 {
937  // Initialise likelihood value
938  double value = 0.0;
939 
940  // Get number of parameters
941  int npars = gradient->size();
942 
943  // Allocate some working arrays
944  int* inx = new int[npars];
945  double* values = new double[npars];
946  GVector wrk_grad(npars);
947 
948  // Determine Npred value and gradient for this observation
949  double npred_value = this->npred(models, &wrk_grad);
950 
951  // Update likelihood, Npred and gradient
952  value += npred_value;
953  *npred += npred_value;
954  *gradient += wrk_grad;
955 
956  // Iterate over all events
957  for (int i = 0; i < events()->size(); ++i) {
958 
959  // Get event pointer
960  const GEvent* event = (*events())[i];
961 
962  // Get model and derivative
963  double model = this->model(models, *event, &wrk_grad);
964 
965  // Skip bin if model is too small (avoids -Inf or NaN gradients)
966  if (model <= minmod) {
967  continue;
968  }
969 
970  // Create index array of non-zero derivatives and initialise working
971  // array
972  int ndev = 0;
973  for (int i = 0; i < npars; ++i) {
974  values[i] = 0.0;
975  if (wrk_grad[i] != 0.0 && !gammalib::is_infinite(wrk_grad[i])) {
976  inx[ndev] = i;
977  ndev++;
978  }
979  }
980 
981  // Update Poissonian statistic (excluding factorial term for faster
982  // computation)
983  value -= log(model);
984 
985  // Skip bin now if there are no non-zero derivatives
986  if (ndev < 1) {
987  continue;
988  }
989 
990  // Update gradient vector and curvature matrix.
991  double fb = 1.0 / model;
992  double fa = fb / model;
993  for (int jdev = 0; jdev < ndev; ++jdev) {
994 
995  // Initialise computation
996  register int jpar = inx[jdev];
997  double g = wrk_grad[jpar];
998  double fa_i = fa * g;
999 
1000  // Update gradient.
1001  (*gradient)[jpar] -= fb * g;
1002 
1003  // Loop over rows
1004  register int* ipar = inx;
1005 
1006  for (register int idev = 0; idev < ndev; ++idev, ++ipar) {
1007  values[idev] = fa_i * wrk_grad[*ipar];
1008  }
1009 
1010  // Add column to matrix
1011  curvature->add_to_column(jpar, values, inx, ndev);
1012 
1013  } // endfor: looped over columns
1014 
1015  } // endfor: iterated over all events
1016 
1017  // Free temporary memory
1018  if (values != NULL) delete [] values;
1019  if (inx != NULL) delete [] inx;
1020 
1021  // Return
1022  return value;
1023 }
1024 
1025 
1026 /***********************************************************************//**
1027  * @brief Evaluate log-likelihood function for Poisson statistic and
1028  * binned analysis (version with working arrays)
1029  *
1030  * @param[in] models Models.
1031  * @param[in,out] gradient Gradient.
1032  * @param[in,out] curvature Curvature matrix.
1033  * @param[in,out] npred Number of predicted events.
1034  * @return Likelihood value.
1035  *
1036  * This method evaluates the -(log-likelihood) function for parameter
1037  * optimisation using binned analysis and Poisson statistic.
1038  * The -(log-likelihood) function is given by
1039  *
1040  * \f[
1041  * L=-\sum_i n_i \log e_i - e_i
1042  * \f]
1043  *
1044  * where the sum is taken over all data space bins, \f$n_i\f$ is the
1045  * observed number of counts and \f$e_i\f$ is the model.
1046  * This method also computes the parameter gradients
1047  * \f$\delta L/dp\f$
1048  * and the curvature matrix
1049  * \f$\delta^2 L/dp_1 dp_2\f$
1050  * and also updates the total number of predicted events m_npred.
1051  ***************************************************************************/
1053  GVector* gradient,
1054  GMatrixSparse* curvature,
1055  double* npred) const
1056 {
1057  // Initialise likelihood value
1058  double value = 0.0;
1059 
1060  // Initialise statistic
1061  #if defined(G_OPT_DEBUG)
1062  int n_bins = 0;
1063  int n_used = 0;
1064  int n_small_model = 0;
1065  int n_zero_data = 0;
1066  int n_exclude_data = 0;
1067  double sum_data = 0.0;
1068  double sum_model = 0.0;
1069  double init_value = value;
1070  #endif
1071 
1072  // Get number of parameters
1073  int npars = gradient->size();
1074 
1075  // Allocate some working arrays
1076  int* inx = new int[npars];
1077  double* values = new double[npars];
1078  GVector wrk_grad(npars);
1079 
1080  // Iterate over all bins
1081  for (int i = 0; i < events()->size(); ++i) {
1082 
1083  // Update number of bins
1084  #if defined(G_OPT_DEBUG)
1085  n_bins++;
1086  #endif
1087 
1088  // Get event pointer
1089  const GEventBin* bin =
1090  (*(static_cast<GEventCube*>(const_cast<GEvents*>(events()))))[i];
1091 
1092  // Get number of counts in bin
1093  double data = bin->counts();
1094 
1095  // Skip bin if data is negative (filtering flag)
1096  if (data < 0) {
1097  #if defined(G_OPT_DEBUG)
1098  n_exclude_data++;
1099  #endif
1100  continue;
1101  }
1102 
1103  // Get model and derivative
1104  double model = this->model(models, *bin, &wrk_grad);
1105 
1106  // Multiply model by bin size
1107  model *= bin->size();
1108 
1109  // Skip bin if model is too small (avoids -Inf or NaN gradients)
1110  if (model <= minmod) {
1111  #if defined(G_OPT_DEBUG)
1112  n_small_model++;
1113  #endif
1114  continue;
1115  }
1116 
1117  // Update statistic
1118  #if defined(G_OPT_DEBUG)
1119  n_used++;
1120  sum_data += data;
1121  sum_model += model;
1122  #endif
1123 
1124  // Update Npred
1125  *npred += model;
1126 
1127  // Multiply gradient by bin size
1128  wrk_grad *= bin->size();
1129 
1130  // Create index array of non-zero derivatives and initialise working
1131  // array
1132  int ndev = 0;
1133  for (int i = 0; i < npars; ++i) {
1134  values[i] = 0.0;
1135  if (wrk_grad[i] != 0.0 && !gammalib::is_infinite(wrk_grad[i])) {
1136  inx[ndev] = i;
1137  ndev++;
1138  }
1139  }
1140 
1141  // Update gradient vector and curvature matrix. To avoid
1142  // unneccessary computations we distinguish the case where
1143  // data>0 and data=0. The second case requires much less
1144  // computation since it does not contribute to the curvature
1145  // matrix ...
1146  if (data > 0.0) {
1147 
1148  // Update Poissonian statistic (excluding factorial term for
1149  // faster computation)
1150  value -= data * log(model) - model;
1151 
1152  // Skip bin now if there are no non-zero derivatives
1153  if (ndev < 1) {
1154  continue;
1155  }
1156 
1157  // Pre computation
1158  double fb = data / model;
1159  double fc = (1.0 - fb);
1160  double fa = fb / model;
1161 
1162  // Loop over columns
1163  for (int jdev = 0; jdev < ndev; ++jdev) {
1164 
1165  // Initialise computation
1166  register int jpar = inx[jdev];
1167  double g = wrk_grad[jpar];
1168  double fa_i = fa * g;
1169 
1170  // Update gradient
1171  (*gradient)[jpar] += fc * g;
1172 
1173  // Loop over rows
1174  register int* ipar = inx;
1175  for (register int idev = 0; idev < ndev; ++idev, ++ipar) {
1176  values[idev] = fa_i * wrk_grad[*ipar];
1177  }
1178 
1179  // Add column to matrix
1180  curvature->add_to_column(jpar, values, inx, ndev);
1181 
1182  } // endfor: looped over columns
1183 
1184  } // endif: data was > 0
1185 
1186  // ... handle now data=0
1187  else {
1188 
1189  // Update statistic
1190  #if defined(G_OPT_DEBUG)
1191  n_zero_data++;
1192  #endif
1193 
1194  // Update Poissonian statistic (excluding factorial term for
1195  // faster computation)
1196  value += model;
1197 
1198  // Skip bin now if there are no non-zero derivatives
1199  if (ndev < 1) {
1200  continue;
1201  }
1202 
1203  // Update gradient
1204  register int* ipar = inx;
1205  for (register int idev = 0; idev < ndev; ++idev, ++ipar) {
1206  (*gradient)[*ipar] += wrk_grad[*ipar];
1207  }
1208 
1209  } // endif: data was 0
1210 
1211  } // endfor: iterated over all events
1212 
1213  // Free temporary memory
1214  if (values != NULL) delete [] values;
1215  if (inx != NULL) delete [] inx;
1216 
1217  // Dump statistic
1218  #if defined(G_OPT_DEBUG)
1219  std::cout << "Number of bins: " << n_bins << std::endl;
1220  std::cout << "Number of bins used for computation: " << n_used << std::endl;
1221  std::cout << "Number of bins excluded from data: " << n_exclude_data << std::endl;
1222  std::cout << "Number of bins excluded due to small model: " << n_small_model << std::endl;
1223  std::cout << "Number of bins with zero data: " << n_zero_data << std::endl;
1224  std::cout << "Sum of data: " << sum_data << std::endl;
1225  std::cout << "Sum of model: " << sum_model << std::endl;
1226  std::cout << "Initial statistic: " << init_value << std::endl;
1227  std::cout << "Statistic: " << value-init_value << std::endl;
1228  #endif
1229 
1230  // Return
1231  return value;
1232 }
1233 
1234 
1235 /***********************************************************************//**
1236  * @brief Evaluate log-likelihood function for Gaussian statistic and
1237  * binned analysis (version with working arrays)
1238  *
1239  * @param[in] models Models.
1240  * @param[in,out] gradient Gradient.
1241  * @param[in,out] curvature Curvature matrix.
1242  * @param[in,out] npred Number of predicted events.
1243  * @return Likelihood value.
1244  *
1245  * This method evaluates the -(log-likelihood) function for parameter
1246  * optimisation using binned analysis and Poisson statistic.
1247  * The -(log-likelihood) function is given by
1248  *
1249  * \f[
1250  * L = 1/2 \sum_i (n_i - e_i)^2 \sigma_i^{-2}
1251  * \f]
1252  *
1253  * where the sum is taken over all data space bins, \f$n_i\f$ is the
1254  * observed number of counts, \f$e_i\f$ is the model and \f$\sigma_i\f$
1255  * is the statistical uncertainty.
1256  * This method also computes the parameter gradients
1257  * \f$\delta L/dp\f$
1258  * and the curvature matrix
1259  * \f$\delta^2 L/dp_1 dp_2\f$
1260  * and also updates the total number of predicted events m_npred.
1261  ***************************************************************************/
1263  GVector* gradient,
1264  GMatrixSparse* curvature,
1265  double* npred) const
1266 {
1267  // Initialise likelihood value
1268  double value = 0.0;
1269 
1270  // Get number of parameters
1271  int npars = gradient->size();
1272 
1273  // Allocate some working arrays
1274  int* inx = new int[npars];
1275  double* values = new double[npars];
1276  GVector wrk_grad(npars);
1277 
1278  // Iterate over all bins
1279  for (int i = 0; i < events()->size(); ++i) {
1280 
1281  // Get event pointer
1282  const GEventBin* bin =
1283  (*(static_cast<GEventCube*>(const_cast<GEvents*>(events()))))[i];
1284 
1285  // Get number of counts in bin
1286  double data = bin->counts();
1287 
1288  // Skip bin if data is negative (filtering flag)
1289  if (data < 0) {
1290  continue;
1291  }
1292 
1293  // Get statistical uncertainty
1294  double sigma = bin->error();
1295 
1296  // Skip bin if statistical uncertainty is too small
1297  if (sigma <= minerr) {
1298  continue;
1299  }
1300 
1301  // Get model and derivative
1302  double model = this->model(models, *bin, &wrk_grad);
1303 
1304  // Multiply model by bin size
1305  model *= bin->size();
1306 
1307  // Skip bin if model is too small (avoids -Inf or NaN gradients)
1308  if (model <= minmod) {
1309  continue;
1310  }
1311 
1312  // Update Npred
1313  *npred += model;
1314 
1315  // Multiply gradient by bin size
1316  wrk_grad *= bin->size();
1317 
1318  // Create index array of non-zero derivatives and initialise working
1319  // array
1320  int ndev = 0;
1321  for (int i = 0; i < npars; ++i) {
1322  values[i] = 0.0;
1323  if (wrk_grad[i] != 0.0 && !gammalib::is_infinite(wrk_grad[i])) {
1324  inx[ndev] = i;
1325  ndev++;
1326  }
1327  }
1328 
1329  // Set weight
1330  double weight = 1.0 / (sigma * sigma);
1331 
1332  // Update Gaussian statistic
1333  double fa = data - model;
1334  value += 0.5 * (fa * fa * weight);
1335 
1336  // Skip bin now if there are no non-zero derivatives
1337  if (ndev < 1) {
1338  continue;
1339  }
1340 
1341  // Loop over columns
1342  for (int jdev = 0; jdev < ndev; ++jdev) {
1343 
1344  // Initialise computation
1345  register int jpar = inx[jdev];
1346  double fa_i = wrk_grad[jpar] * weight;
1347 
1348  // Update gradient
1349  (*gradient)[jpar] -= fa * fa_i;
1350 
1351  // Loop over rows
1352  register int* ipar = inx;
1353  for (register int idev = 0; idev < ndev; ++idev, ++ipar) {
1354  values[idev] = fa_i * wrk_grad[*ipar];
1355  }
1356 
1357  // Add column to matrix
1358  curvature->add_to_column(jpar, values, inx, ndev);
1359 
1360  } // endfor: looped over columns
1361 
1362  } // endfor: iterated over all events
1363 
1364  // Free temporary memory
1365  if (values != NULL) delete [] values;
1366  if (inx != NULL) delete [] inx;
1367 
1368  // Return
1369  return value;
1370 }
1371 
1372 
1373 /*==========================================================================
1374  = =
1375  = Model gradient methods =
1376  = =
1377  ==========================================================================*/
1378 
1379 /***********************************************************************//**
1380  * @brief Model function evaluation for gradient computation
1381  *
1382  * @param[in] x Function value.
1383  ***************************************************************************/
1384 double GObservation::model_func::eval(const double& x)
1385 {
1386  // Set value
1387  m_par->factor_value(x);
1388 
1389  // Compute model value
1390  double value = m_model->eval(*m_event, *m_parent);
1391 
1392  // Return value
1393  return value;
1394 }
1395 
1396 
1397 /*==========================================================================
1398  = =
1399  = Npred computation methods =
1400  = =
1401  ==========================================================================*/
1402 
1403 /***********************************************************************//**
1404  * @brief Npred function evaluation for gradient computation
1405  *
1406  * @param[in] x Function value.
1407  ***************************************************************************/
1408 double GObservation::npred_func::eval(const double& x)
1409 {
1410  // Set value
1411  m_par->factor_value(x);
1412 
1413  // Compute Npred value
1414  double npred = m_parent->npred(*m_model);
1415 
1416  // Return value
1417  return npred;
1418 }
1419 
1420 
1421 /***********************************************************************//**
1422  * @brief Integration kernel for npred() method
1423  *
1424  * @param[in] x Function value.
1425  ***************************************************************************/
1426 double GObservation::npred_kern::eval(const double& x)
1427 {
1428  // Convert argument in native reference in seconds
1429  GTime time;
1430  time.secs(x);
1431 
1432  // Return value
1433  return (m_parent->npred_spec(*m_model, time));
1434 }
1435 
1436 
1437 /***********************************************************************//**
1438  * @brief Integrates spatially integrated Npred kernel spectrally
1439  *
1440  * @param[in] model Gamma-ray source model.
1441  * @param[in] obsTime Measured photon arrival time.
1442  *
1443  * @exception GException::erange_invalid
1444  * Energy range is invalid.
1445  *
1446  * Computes
1447  *
1448  * \f[
1449  * N_{\rm pred}(t') = \int_{E_{\rm bounds}} \int_{\rm ROI}
1450  * P(p',E',t') \, dp' \, dE'
1451  * \f]
1452  *
1453  * of the event probability
1454  *
1455  * \f[
1456  * P(p',E',t') = \int \int \int
1457  * S(p,E,t) \times R(p',E',t'|p,E,t) \, dp \, dE \, dt
1458  * \f]
1459  *
1460  * where
1461  * \f$S(p,E,t)\f$ is the source model,
1462  * \f$R(p',E',t'|p,E,t)\f$ is the instrument response function,
1463  * \f$p'\f$ is the measured photon direction,
1464  * \f$E'\f$ is the measured photon energy,
1465  * \f$t'\f$ is the measured photon arrival time,
1466  * \f$p\f$ is the true photon arrival direction,
1467  * \f$E\f$ is the true photon energy, and
1468  * \f$t\f$ is the true photon arrival time.
1469  *
1470  * This method performs the energy intergration over the energy boundaries
1471  * \f$E_{\rm bounds}\f$.
1472  ***************************************************************************/
1474  const GTime& obsTime) const
1475 {
1476  // Set number of iterations for Romberg integration.
1477  static const int iter = 8;
1478 
1479  // Initialise result
1480  double result = 0.0;
1481 
1482  // Get energy boundaries
1483  GEbounds ebounds = events()->ebounds();
1484 
1485  // Loop over energy boundaries
1486  for (int i = 0; i < ebounds.size(); ++i) {
1487 
1488  // Get boundaries in MeV
1489  double emin = ebounds.emin(i).MeV();
1490  double emax = ebounds.emax(i).MeV();
1491 
1492  // Continue only if valid
1493  if (emax > emin) {
1494 
1495  // Setup integration function
1496  GObservation::npred_spec_kern integrand(this, &model, &obsTime);
1497  GIntegral integral(&integrand);
1498 
1499  // Set number of iterations
1500  integral.fixed_iter(iter);
1501 
1502  // Do Romberg integration
1503  #if defined(G_LN_ENERGY_INT)
1504  emin = std::log(emin);
1505  emax = std::log(emax);
1506  #endif
1507  result += integral.romberg(emin, emax);
1508 
1509  } // endif: energy interval was valid
1510 
1511  } // endfor: looped over energy boundaries
1512 
1513  // Compile option: Check for NaN/Inf
1514  #if defined(G_NAN_CHECK)
1515  if (gammalib::is_notanumber(result) || gammalib::is_infinite(result)) {
1516  std::cout << "*** ERROR: GObservation::npred_spec:";
1517  std::cout << " NaN/Inf encountered";
1518  std::cout << " (result=" << result;
1519  std::cout << ")" << std::endl;
1520  }
1521  #endif
1522 
1523  // Return result
1524  return result;
1525 }
1526 
1527 
1528 /***********************************************************************//**
1529  * @brief Integration kernel for npred_spec() method
1530  *
1531  * @param[in] x Function value.
1532  *
1533  * This method implements the integration kernel needed for the npred_spec()
1534  * method. If G_LN_ENERGY_INT is defined the energy integration is done
1535  * logarithmically, i.e. @p x is given in ln(energy) instead of energy.
1536  ***************************************************************************/
1538 {
1539  // Set energy
1540  GEnergy eng;
1541  #if defined(G_LN_ENERGY_INT)
1542  double expx = std::exp(x);
1543  eng.MeV(expx);
1544  #else
1545  eng.MeV(x);
1546  #endif
1547 
1548  // Get function value
1549  double value = m_model->npred(eng, *m_time, *m_parent);
1550 
1551  // Save value if needed
1552  #if defined(G_NAN_CHECK)
1553  double value_out = value;
1554  #endif
1555 
1556  // Correct for variable substitution
1557  #if defined(G_LN_ENERGY_INT)
1558  value *= expx;
1559  #endif
1560 
1561  // Compile option: Check for NaN
1562  #if defined(G_NAN_CHECK)
1563  if (gammalib::is_notanumber(value) || gammalib::is_infinite(value)) {
1564  std::cout << "*** ERROR: GObservation::npred_spec_kern::eval";
1565  std::cout << "(x=" << x << "): ";
1566  std::cout << " NaN/Inf encountered";
1567  std::cout << " (value=" << value;
1568  std::cout << " (value_out=" << value_out;
1569  #if defined(G_LN_ENERGY_INT)
1570  std::cout << " exp(x)=" << expx;
1571  #endif
1572  std::cout << ")" << std::endl;
1573  }
1574  #endif
1575 
1576  // Return value
1577  return value;
1578 }
virtual double error(void) const =0
Return error in number of counts.
Definition: GEventBin.cpp:141
const std::string & statistic(void) const
Return optimizer statistic.
Abstract model class.
Definition: GModel.hpp:97
double eval(const double &x)
Integration kernel for npred_spec() method.
virtual double likelihood_gaussian_binned(const GModels &models, GVector *gradient, GMatrixSparse *curvature, double *npred) const
Evaluate log-likelihood function for Gaussian statistic and binned analysis (version with working arr...
virtual bool is_constant(void) const =0
const GEvent * m_event
Event.
virtual double likelihood_poisson_unbinned(const GModels &models, GVector *gradient, GMatrixSparse *curvature, double *npred) const
Evaluate log-likelihood function for Poisson statistic and unbinned analysis (version with working ar...
Abstract model base class interface definition.
double eval(const double &x)
Integration kernel for npred() method.
virtual double size(void) const =0
const double & factor_gradient(void) const
Return parameter gradient factor.
virtual double ontime(void) const =0
double romberg(std::vector< double > bounds, const int &order=5)
Perform Romberg integration.
Definition: GIntegral.cpp:380
std::string m_name
Observation name.
Sparse matrix class interface definition.
virtual std::string instrument(void) const =0
virtual double npred(const GModels &models, GVector *gradient=NULL) const
Return total number (and optionally gradient) of predicted counts for all models. ...
Model container class definition.
void ebounds(const GEbounds &ebounds)
Set energy boundaries.
Definition: GEvents.cpp:136
virtual double likelihood(const GModels &models, GVector *gradient, GMatrixSparse *curvature, double *npred) const
Compute likelihood function.
int size(void) const
Return number of energy boundaries.
Definition: GEbounds.hpp:157
virtual double npred_grad(const GModel &model, const GModelPar &par) const
Returns parameter gradient of Npred.
GEvents * m_events
Pointer to event container.
Abstract interface for the event bin class.
Definition: GEventBin.hpp:64
virtual double eval(const GEvent &event, const GObservation &obs, const bool &gradients=false) const =0
Abstract event bin base class definition.
Abstract event bin container class interface definition.
Abstract interface for the event classes.
Definition: GEvent.hpp:71
void gti(const GGti &gti)
Set Good Time Intervals.
Definition: GEvents.cpp:154
virtual double counts(void) const =0
std::string m_id
Observation identifier.
double MeV(void) const
Return energy in MeV.
Definition: GEnergy.cpp:321
virtual double model_grad(const GModel &model, const GModelPar &par, const GEvent &event) const
Returns parameter gradient of model for a given event.
Time class.
Definition: GTime.hpp:54
Gammalib tools definition.
GIntegral class interface definition.
Definition: GIntegral.hpp:45
virtual void remove_response_cache(const std::string &name)
Response cache removal hook.
virtual double npred_spec(const GModel &model, const GTime &obsTime) const
Integrates spatially integrated Npred kernel spectrally.
bool is_free(void) const
Signal if parameter is free.
virtual int nobserved(void) const
Return total number of observed events.
Abstract event atom container class interface definition.
void free_members(void)
Delete class members.
virtual void add_to_column(const int &column, const GVector &vector)
Add vector column into matrix.
virtual ~GObservation(void)
Destructor.
bool is_notanumber(const double &x)
Signal if argument is not a number.
Definition: GTools.hpp:178
Model parameter class interface definition.
int size(void) const
Return number of parameters in model.
Definition: GModel.hpp:217
bool is_infinite(const double &x)
Signal if argument is infinite.
Definition: GTools.hpp:161
const double minerr
Minimum statistical error.
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
Model parameter class.
Definition: GModelPar.hpp:87
#define G_MODEL
Model container class.
Definition: GModels.hpp:150
virtual double model(const GModels &models, const GEvent &event, GVector *gradient=NULL) const
Return model value and (optionally) gradient.
virtual double likelihood_poisson_binned(const GModels &models, GVector *gradient, GMatrixSparse *curvature, double *npred) const
Evaluate log-likelihood function for Poisson statistic and binned analysis (version with working arra...
bool has_grad(void) const
Signal if parameter gradient is computed analytically.
#define G_NPRED
Energy boundaries container class.
Definition: GEbounds.hpp:60
const double & factor_max(void) const
Return parameter maximum boundary factor.
const GEnergy & emin(void) const
Return minimum energy of all intervals.
Definition: GEbounds.hpp:181
Abstract data model base class interface definition.
GObservation(void)
Void constructor.
#define G_EVENTS
GVector log(const GVector &vector)
Computes natural logarithm of vector elements.
Definition: GVector.cpp:1184
double eval(const double &x)
Model function evaluation for gradient computation.
GModelPar * m_par
Model parameter.
bool has_max(void) const
Signal if parameter has maximum boundary.
void init_members(void)
Initialise class members.
Abstract observation base class.
Vector class interface definition.
Abstract observation base class interface definition.
const double minmod
Minimum model value.
virtual GEvents * clone(void) const =0
Clones object.
virtual bool use_edisp(void) const =0
Abstract response base class definition.
const double & factor_min(void) const
Return parameter minimum boundary factor.
const double & secs(void) const
Return time in seconds in native reference (TT)
Definition: GTime.hpp:153
#define G_LIKELIHOOD
Sky model class interface definition.
Abstract event container class.
Definition: GEvents.hpp:66
Numerical derivatives class.
Definition: GDerivative.hpp:45
GDerivative class interface definition.
Sparse matrix class definition.
bool has_min(void) const
Signal if parameter has minimum boundary.
std::string m_statistic
Optimizer statistic.
virtual GObservation & operator=(const GObservation &obs)
Assignment operator.
virtual GEvents * events(void)
Return events.
Exception handler interface definition.
int npars(void) const
Return number of model parameters in container.
Definition: GModels.cpp:813
const GModel * m_model
Model.
void fixed_iter(const int &iter)
Set fixed number of iterations.
Definition: GIntegral.hpp:181
std::string toupper(const std::string &s)
Convert string to upper case.
Definition: GTools.cpp:818
GVector exp(const GVector &vector)
Computes exponential of vector elements.
Definition: GVector.cpp:1142
const int & size(void) const
Return size of vector.
Definition: GVector.hpp:180
double eval(const double &x)
Npred function evaluation for gradient computation.
Vector class.
Definition: GVector.hpp:46
const GEnergy & emax(void) const
Return maximum energy of all intervals.
Definition: GEbounds.hpp:193
void copy_members(const GObservation &obs)
Copy class members.
virtual int number(void) const =0
Integration class interface definition.
virtual const GResponse * response(void) const =0
Abstract event bin container class.
Definition: GEventCube.hpp:46
int size(void) const
Return number of models in container.
Definition: GModels.hpp:255
const GObservation * m_parent
Observation.
const double & factor_value(void) const
Return parameter value factor.
Class that handles energies in a unit independent way.
Definition: GEnergy.hpp:48
virtual int size(void) const =0
std::string str(const unsigned short int &value)
Convert unsigned short integer value into string.
Definition: GTools.cpp:411