GammaLib  2.0.0
 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-2021 by Juergen Knoedlseder *
5  * ----------------------------------------------------------------------- *
6  * *
7  * This program is free software: you can redistribute it and/or modify *
8  * it under the terms of the GNU General Public License as published by *
9  * the Free Software Foundation, either version 3 of the License, or *
10  * (at your option) any later version. *
11  * *
12  * This program is distributed in the hope that it will be useful, *
13  * but WITHOUT ANY WARRANTY; without even the implied warranty of *
14  * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the *
15  * GNU General Public License for more details. *
16  * *
17  * You should have received a copy of the GNU General Public License *
18  * along with this program. If not, see <http://www.gnu.org/licenses/>. *
19  * *
20  ***************************************************************************/
21 /**
22  * @file 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_MODEL1 "GObservation::model(GModels&, GEvent&, GVector*)"
52 #define G_MODEL2 "GObservation::model(GModels&, GMatrixSparse*)"
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 //!< Debug likelihood computation
68 //#define G_DEBUG_VECTOR_MODEL //!< Debug vector model
69 
70 
71 /*==========================================================================
72  = =
73  = Constructors/destructors =
74  = =
75  ==========================================================================*/
76 
77 /***********************************************************************//**
78  * @brief Void constructor
79  ***************************************************************************/
81 {
82  // Initialise members
83  init_members();
84 
85  // Return
86  return;
87 }
88 
89 
90 /***********************************************************************//**
91  * @brief Copy constructor
92  *
93  * @param[in] obs Observation.
94  *
95  * Instantiate the class by copying from an existing observation.
96  ***************************************************************************/
98 {
99  // Initialise members
100  init_members();
101 
102  // Copy members
103  copy_members(obs);
104 
105  // Return
106  return;
107 }
108 
109 
110 /***********************************************************************//**
111  * @brief Destructor
112  ***************************************************************************/
114 {
115  // Free members
116  free_members();
117 
118  // Return
119  return;
120 }
121 
122 
123 /*==========================================================================
124  = =
125  = Operators =
126  = =
127  ==========================================================================*/
128 
129 /***********************************************************************//**
130  * @brief Assignment operator
131  *
132  * @param[in] obs Observation.
133  * @return Observation.
134  *
135  * Assign observation.
136  ***************************************************************************/
138 {
139  // Execute only if object is not identical
140  if (this != &obs) {
141 
142  // Free members
143  free_members();
144 
145  // Initialise members
146  init_members();
147 
148  // Copy members
149  copy_members(obs);
150 
151  } // endif: object was not identical
152 
153  // Return this object
154  return *this;
155 }
156 
157 
158 /*==========================================================================
159  = =
160  = Public methods =
161  = =
162  ==========================================================================*/
163 
164 /***********************************************************************//**
165  * @brief Return events
166  *
167  * @exception GException::no_events
168  * No events allocated for observation.
169  *
170  * Returns pointer to events.
171  ***************************************************************************/
173 {
174  // Throw an exception if the event container is not valid
175  if (m_events == NULL) {
176  std::string msg = "No events allocated for observation.";
178  }
179 
180  // Return pointer to events
181  return (m_events);
182 }
183 
184 
185 /***********************************************************************//**
186  * @brief Return events (const version)
187  *
188  * @exception GException::no_events
189  * No events allocated for observation.
190  *
191  * Returns const pointer to events.
192  ***************************************************************************/
193 const GEvents* GObservation::events(void) const
194 {
195  // Throw an exception if the event container is not valid
196  if (m_events == NULL) {
197  std::string msg = "No events allocated for observation.";
199  }
200 
201  // Return pointer to events
202  return (m_events);
203 }
204 
205 
206 /***********************************************************************//**
207  * @brief Set event container
208  *
209  * @param[in] events Event container.
210  *
211  * Set the event container for this observation by cloning the container
212  * specified in the argument.
213  ***************************************************************************/
214 void GObservation::events(const GEvents& events)
215 {
216  // Free existing event container only if it differs from current event
217  // container. This prevents unintential deallocation of the argument
218  if ((m_events != NULL) && (m_events != &events)) {
219  delete m_events;
220  }
221 
222  // Clone events
223  m_events = events.clone();
224 
225  // Return
226  return;
227 }
228 
229 
230 /***********************************************************************//**
231  * @brief Compute likelihood function
232  *
233  * @param[in] models Models.
234  * @param[in,out] gradients Pointer to gradients.
235  * @param[in,out] curvature Pointer to curvature matrix.
236  * @param[in,out] npred Pointer to Npred value.
237  * @return Likelihood.
238  *
239  * Computes the likelihood for a specified set of models. The method also
240  * returns the gradients, the curvature matrix, and the number of events
241  * that are predicted by all models.
242  ***************************************************************************/
243 double GObservation::likelihood(const GModels& models,
244  GVector* gradients,
245  GMatrixSparse* curvature,
246  double* npred) const
247 {
248  // Initialise likelihood value
249  double value = 0.0;
250 
251  // Extract statistic for this observation
252  std::string statistic = gammalib::toupper(this->statistic());
253 
254  // Unbinned analysis
255  if (dynamic_cast<const GEventList*>(events()) != NULL) {
256 
257  // Poisson statistic
258  if ((statistic == "POISSON") || (statistic == "CSTAT")) {
259 
260  // Update the log-likelihood
261  value = likelihood_poisson_unbinned(models,
262  gradients,
263  curvature,
264  npred);
265 
266  } // endif: Poisson statistic
267 
268  // ... otherwise throw an exception
269  else {
270  std::string msg = "Invalid statistic \""+statistic+"\". Unbinned "
271  "optimization requires \"POISSON\" or \"CSTAT\" "
272  "statistic.";
274  }
275 
276  } // endif: unbinned analysis
277 
278  // ... or binned analysis
279  else {
280 
281  // Poisson statistic
282  if ((statistic == "POISSON") || (statistic == "CSTAT")) {
283  value = likelihood_poisson_binned(models,
284  gradients,
285  curvature,
286  npred);
287  }
288 
289  // ... or Gaussian statistic
290  else if ((statistic == "GAUSSIAN") || (statistic == "CHI2")) {
291  value = likelihood_gaussian_binned(models,
292  gradients,
293  curvature,
294  npred);
295  }
296 
297  // ... or unsupported
298  else {
299  std::string msg = "Invalid statistic \""+statistic+"\". Binned "
300  "optimization requires \"POISSON\", \"CSTAT\", "
301  "\"GAUSSIAN\" or \"CHI2\" statistic.";
303  }
304 
305  } // endelse: binned analysis
306 
307  // Return likelihood
308  return value;
309 }
310 
311 
312 /***********************************************************************//**
313  * @brief Return model value and (optionally) gradients
314  *
315  * @param[in] models Model container.
316  * @param[in] event Observed event.
317  * @param[out] gradients Pointer to gradient vector (optional).
318  * @return Model value.
319  *
320  * @exception GException::invalid_value
321  * Dimension of gradient vector mismatches number of parameters.
322  *
323  * Implements generic model and gradient evaluation for an observation. The
324  * model gives the probability for an event to occur with a given instrument
325  * direction, at a given energy and at a given time. The gradient is the
326  * parameter derivative of this probability. If NULL is passed for the
327  * gradient vector, then gradients will not be computed.
328  *
329  * The method will only operate on models for which the list of instruments
330  * and observation identifiers matches those of the observation. Models that
331  * do not match will be skipped.
332  ***************************************************************************/
333 double GObservation::model(const GModels& models,
334  const GEvent& event,
335  GVector* gradients) const
336 {
337  // Initialise method variables
338  double model = 0.0; // Reset model value
339  int igrad = 0; // Reset gradient counter
340  int grad_size = 0; // Reset gradient size
341 
342  // Set gradient usage flag
343  bool use_grad = (gradients != NULL);
344 
345  // If gradient is available then reset gradient vector elements to 0
346  // and determine vector size
347  if (use_grad) {
348 
349  // Reset gradient vector
350  (*gradients) = 0.0;
351  grad_size = gradients->size();
352 
353  // Initialise stack of parameters with gradients
354  m_pars_with_gradients.clear();
355  m_pars_with_gradients.reserve(models.size());
356 
357  } // endif: gradient vector available
358 
359  // Loop over models
360  for (int i = 0; i < models.size(); ++i) {
361 
362  // Get model pointer. Continue only if pointer is valid
363  const GModel* mptr = models[i];
364  if (mptr != NULL) {
365 
366  // Continue only if model applies to specific instrument and
367  // observation identifier
368  if (mptr->is_valid(instrument(), id())) {
369 
370  // Compute value and add to model
371  model += mptr->eval(event, *this, use_grad);
372 
373  // Optionally determine model gradients
374  if (use_grad) {
375 
376  // Make sure that we have a slot for the gradient
377  #if defined(G_RANGE_CHECK)
378  if (igrad+mptr->size() > grad_size) {
379  std::string msg = "Vector has not enough elements to "
380  "store the model parameter gradients. "+
381  gammalib::str(models.npars())+" elements "
382  "requested while vector only contains "+
383  gammalib::str(gradients->size())+" elements.";
385  }
386  #endif
387 
388  // If model provides the parameter indices that were
389  // updated by the eval() method then use them as this is
390  // less time consuming than looping over all indices
391  if (mptr->has_eval_indices()) {
392 
393  // Get number of relevant parameters
394  int npars = mptr->eval_indices().size();
395 
396  // Loop over all relevant parameters
397  for (int index = 0; index < npars; ++index) {
398 
399  // Retrieve parameter index
400  int ipar = mptr->eval_indices()[index];
401 
402  // Get reference to model parameter
403  const GModelPar& par = (*mptr)[ipar];
404 
405  // Set gradient
406  if (par.is_free()) {
407  if (has_gradient(*mptr, par)) {
408  (*gradients)[igrad+ipar] =
409  par.factor_gradient();
410  }
411  else {
412  (*gradients)[igrad+ipar] =
413  model_grad(*mptr, par, event);
414  }
415  }
416  else {
417  (*gradients)[igrad+ipar] = 0.0;
418  }
419 
420  } // endfor: looped over all parameter indices
421 
422  } // endif: parameter indices were available
423 
424  // ... otherwise loop over all parameter indices
425  else {
426 
427  // Loop over all parameters
428  for (int ipar = 0; ipar < mptr->size(); ++ipar) {
429 
430  // Get reference to model parameter
431  const GModelPar& par = (*mptr)[ipar];
432 
433  // Set gradient
434  if (par.is_free()) {
435  if (has_gradient(*mptr, par)) {
436  (*gradients)[igrad+ipar] =
437  par.factor_gradient();
438  }
439  else {
440  (*gradients)[igrad+ipar] =
441  model_grad(*mptr, par, event);
442  }
443  }
444  else {
445  (*gradients)[igrad+ipar] = 0.0;
446  }
447 
448  } // endfor: looped over all parameter indices
449 
450  } // endelse: no parameter indices were available
451 
452  } // endif: use model gradients
453 
454  } // endif: model component was valid for instrument
455 
456  // Increment parameter counter for gradients
457  igrad += mptr->size();
458 
459  } // endif: model was valid
460 
461  } // endfor: Looped over models
462 
463  // Return
464  return model;
465 }
466 
467 
468 /***********************************************************************//**
469  * @brief Return vector of model values and (optionally) gradients
470  *
471  * @param[in] models Model container.
472  * @param[out] gradients Pointer to sparse gradient matrix.
473  * @return Vector of model values.
474  *
475  * @exception GException::invalid_argument
476  * Gradient matrix mismatches number of events or model parameters.
477  *
478  * Returns the model values for each event in the observation.
479  *
480  * If @p gradients is not NULL, the matrix contains on output the model
481  * factor gradients for all events. Each row of the @p gradients matrix
482  * corresponds to one event, the columns correspond to the parameters
483  * of the @p models container.
484  ***************************************************************************/
486  GMatrixSparse* gradients) const
487 {
488  // Initialise variables
489  int nevents = events()->size();
490  bool use_grad = (gradients != NULL);
491  int igrad = 0;
492 
493  // Initialise model values
494  GVector values(nevents);
495 
496  // If gradient is available then check gradient size and initialise
497  // sparse matrix stack
498  if (use_grad) {
499 
500  // Initialise stack of parameters with gradients
501  m_pars_with_gradients.clear();
502  m_pars_with_gradients.reserve(models.size());
503 
504  // Check number of columns
505  int ncolumns = (nevents > 0) ? models.npars() : 0;
506  if (gradients->columns() != ncolumns) {
507  std::string msg = "Number of "+gammalib::str(gradients->columns())+" "
508  "columns in gradient matrix differs from expected "
509  "number of "+gammalib::str(ncolumns)+". Please "
510  "specify compatible arguments.";
512  }
513 
514  // Check number of rows
515  int nrows = (ncolumns > 0) ? nevents : 0;
516  if (gradients->rows() != nrows) {
517  std::string msg = "Number of "+gammalib::str(gradients->rows())+" "
518  "rows in gradient matrix differs from expected "
519  "number of "+gammalib::str(nrows)+". Please "
520  "specify compatible arguments.";
522  }
523 
524  } // endif: gradient was requested
525 
526  // Loop over models
527  for (int i = 0; i < models.size(); ++i) {
528 
529  // Get model pointer. Continue only if pointer is valid
530  const GModel* mptr = models[i];
531  if (mptr != NULL) {
532 
533  // Continue only if model applies to specific instrument and
534  // observation identifier
535  if (mptr->is_valid(instrument(), id())) {
536 
537  // If no gradients should be used then evaluate model and
538  // add it to values
539  if (!use_grad) {
540  values += mptr->eval(*this, NULL);
541  }
542 
543  // ... otherwise if there are events then evaluate model and
544  // gradients and add to model vector and gradient matrix
545  else if (nevents > 0) {
546 
547  // Allocate gradient matrix
548  GMatrixSparse grads(nevents, mptr->size());
549 
550  // Initialise sparse matrix stack
551  grads.stack_init(nevents*mptr->size(), mptr->size());
552 
553  // Evaluate model and add to values
554  values += mptr->eval(*this, &grads);
555 
556  // Destroy sparse matrix stack (fills stack in matrix)
557  grads.stack_destroy();
558 
559  // If model provides the parameter indices that were
560  // updated by the eval() method then use them as this
561  // is less time consuming than looping over all indices
562  if (mptr->has_eval_indices()) {
563 
564  // Get number of relevant parameters
565  int npars = mptr->eval_indices().size();
566 
567  // Loop over all relevant parameters
568  for (int index = 0; index < npars; ++index) {
569 
570  // Retrieve parameter index
571  int ipar = mptr->eval_indices()[index];
572 
573  // Get reference to model parameter
574  const GModelPar& par = (*mptr)[ipar];
575 
576  // If parameters is free then set gradient
577  if (par.is_free()) {
578 
579  // Determine gradient
580  GVector grad(nevents);
581  if (has_gradient(*mptr, par)) {
582  grad = grads.column(ipar);
583  }
584  else {
585  grad = model_grad(*mptr, par);
586  }
587 
588  // Set gradient
589  gradients->column(igrad+ipar, grad);
590 
591  } // endif: parameter was free
592 
593  } // endfor: looped over all parameter indices
594 
595  } // endif: parameter indices were available
596 
597  // ... otherwise loop over all parameter indices
598  else {
599 
600  // Loop over all parameters
601  for (int ipar = 0; ipar < mptr->size(); ++ipar) {
602 
603  // Get reference to model parameter
604  const GModelPar& par = (*mptr)[ipar];
605 
606  // If parameters is free then set gradient
607  if (par.is_free()) {
608 
609  // Determine gradient
610  GVector grad(nevents);
611  if (has_gradient(*mptr, par)) {
612  grad = grads.column(ipar);
613  #if defined(G_DEBUG_VECTOR_MODEL)
614  if (ipar < 3) {
615  GVector num_grad = model_grad(*mptr, par);
616  for (int k = 0; k < nevents; ++k) {
617  if (grad[k] != 0.0 || num_grad[k] != 0.0) {
618  std::cout << "ipar=" << ipar;
619  std::cout << " par=" << par.name();
620  std::cout << " k=" << k;
621  std::cout << " ana=" << grad[k];
622  std::cout << " num=" << num_grad[k];
623  std::cout << std::endl;
624  }
625  }
626  }
627  #endif
628  }
629  else {
630  grad = model_grad(*mptr, par);
631  }
632 
633  // Set gradient
634  gradients->column(igrad+ipar, grad);
635 
636  } // endif: parameter was free
637 
638  } // endfor: looped over all parameter indices
639 
640  } // endelse: no parameter indices were available
641 
642  } // endif: use model gradients
643 
644  } // endif: model component was valid for instrument
645 
646  // Increment parameter counter for gradients
647  igrad += mptr->size();
648 
649  } // endif: model was valid
650 
651  } // endfor: Looped over models
652 
653  // Return values
654  return values;
655 }
656 
657 
658 /***********************************************************************//**
659  * @brief Return total number of observed events
660  *
661  * @returns Total number of observed events.
662  *
663  * Returns the total number of observed events. If the observation does not
664  * contain any events the method returns zero.
665  ***************************************************************************/
666 int GObservation::nobserved(void) const
667 {
668  // Initialise number of observed events
669  int nobserved = 0;
670 
671  // Extract number of observed events
672  if (m_events != NULL) {
673  nobserved = m_events->number();
674  }
675 
676  // Return number of observed events
677  return nobserved;
678 }
679 
680 
681 /***********************************************************************//**
682  * @brief Return total number (and optionally gradients) of predicted counts
683  * for all models
684  *
685  * @param[in] models Models.
686  * @param[out] gradients Model parameter gradients (optional).
687  *
688  * @exception GException::invalid_argument
689  * Dimension of gradient vector mismatches number of parameters.
690  *
691  * Returns the total number of predicted counts within the analysis region.
692  * If NULL is passed for the gradient vector then gradients will not be
693  * computed.
694  *
695  * The method will only operate on models for which the list of instruments
696  * and observation identifiers matches those of the observation. Models that
697  * do not match will be skipped.
698  ***************************************************************************/
699 double GObservation::npred(const GModels& models, GVector* gradients) const
700 {
701  // Verify that gradient vector and models have the same dimension
702  #if defined(G_RANGE_CHECK)
703  if (gradients != NULL) {
704  if (models.npars() != gradients->size()) {
705  std::string msg = "Model parameter gradients vector size "+
706  gammalib::str(gradients->size())+" differs from "
707  "number of "+gammalib::str(models.npars())+
708  " model parameters. Please specify a gradients "
709  "vector with the appropriate size.";
711  }
712  }
713  #endif
714 
715  // Initialise
716  double npred = 0.0; // Reset predicted number of counts
717  int igrad = 0; // Reset gradient counter
718 
719  // If gradient is available then reset gradient vector elements to 0
720  if (gradients != NULL) {
721  (*gradients) = 0.0;
722  }
723 
724  // Loop over models
725  for (int i = 0; i < models.size(); ++i) {
726 
727  // Get model pointer. Continue only if pointer is valid
728  const GModel* mptr = models[i];
729  if (mptr != NULL) {
730 
731  // Continue only if model applies to specific instrument and
732  // observation identifier
733  if (mptr->is_valid(instrument(), id())) {
734 
735  // Determine Npred for model
736  npred += this->npred(*mptr);
737 
738  // Optionally determine Npred gradients
739  if (gradients != NULL) {
740  for (int k = 0; k < mptr->size(); ++k) {
741  const GModelPar& par = (*mptr)[k];
742  (*gradients)[igrad+k] = npred_grad(*mptr, par);
743  }
744  }
745 
746  } // endif: model component was valid for instrument
747 
748  // Increment parameter counter for gradient
749  igrad += mptr->size();
750 
751  } // endif: model was valid
752 
753  } // endfor: Looped over models
754 
755  // Return prediction
756  return npred;
757 }
758 
759 
760 /***********************************************************************//**
761  * @brief Return total number of predicted counts for one model
762  *
763  * @param[in] model Gamma-ray source model.
764  *
765  * @exception GException::gti_invalid
766  * Good Time Interval is invalid.
767  *
768  * Computes
769  *
770  * \f[
771  * N_{\rm pred} = \int_{\rm GTI} \int_{E_{\rm bounds}} \int_{\rm ROI}
772  * P(p',E',t') \, dp' \, dE' \, dt'
773  * \f]
774  *
775  * of the event probability
776  *
777  * \f[
778  * P(p',E',t') = \int \int \int
779  * S(p,E,t) \times R(p',E',t'|p,E,t) \, dp \, dE \, dt
780  * \f]
781  *
782  * where
783  * \f$S(p,E,t)\f$ is the source model,
784  * \f$R(p',E',t'|p,E,t)\f$ is the instrument response function,
785  * \f$p'\f$ is the measured photon direction,
786  * \f$E'\f$ is the measured photon energy,
787  * \f$t'\f$ is the measured photon arrival time,
788  * \f$p\f$ is the true photon arrival direction,
789  * \f$E\f$ is the true photon energy, and
790  * \f$t\f$ is the true photon arrival time.
791  *
792  * This method performs the time integration over the Good Time Intervals
793  * \f${\rm GTI}\f$. Note that MET is used for the time integration interval.
794  * This, however, is no specialisation since npred_grad_kern_spec::eval()
795  * converts the argument back in a GTime object by assuming that the argument
796  * is in MET, hence the correct time system will be used at the end by the
797  * method.
798  ***************************************************************************/
799 double GObservation::npred(const GModel& model) const
800 {
801  // Initialise result
802  double npred = 0.0;
803 
804  // Continue only if model applies to specific instrument and
805  // observation identifier
806  if (model.is_valid(instrument(), id())) {
807 
808  // Case A: If the model is constant then integrate analytically
809  if (model.is_constant()) {
810 
811  // Evaluate model at first start time and multiply by ontime
812  double ontime = events()->gti().ontime();
813 
814  // Integrate only if ontime is positive
815  if (ontime > 0.0) {
816 
817  // Integration is a simple multiplication by the time
818  npred = npred_spec(model, events()->gti().tstart()) * ontime;
819 
820  }
821 
822  } // endif: model was constant
823 
824  // ... otherwise integrate temporally
825  else {
826 
827  // Loop over GTIs
828  for (int i = 0; i < events()->gti().size(); ++i) {
829 
830  // Set integration interval in seconds
831  double tstart = events()->gti().tstart(i).secs();
832  double tstop = events()->gti().tstop(i).secs();
833 
834  // Throw exception if time interval is not valid
835  if (tstop <= tstart) {
836  std::string msg = "Start time "+
837  events()->gti().tstart(i).print()+
838  " is past the stop time "+
839  events()->gti().tstop(i).print()+
840  " for Good Time Interval "+
841  gammalib::str(i)+". Please make sure "
842  "that all Good Time Intervals have start "
843  "times that are earlier than the stop "
844  "time.";
846  }
847 
848  // Setup integration function
849  GObservation::npred_kern integrand(this, &model);
850  GIntegral integral(&integrand);
851 
852  // Do Romberg integration
853  npred += integral.romberg(tstart, tstop);
854 
855  } // endfor: looped over GTIs
856 
857  } // endelse: integrated temporally
858 
859  } // endif: model was valid
860 
861  // Return Npred
862  return npred;
863 }
864 
865 
866 /***********************************************************************//**
867  * @brief Returns parameter gradient of model for a given event
868  *
869  * @param[in] model Model.
870  * @param[in] par Model parameter.
871  * @param[in] event Event.
872  *
873  * This method uses a robust but simple difference method to estimate
874  * parameter gradients that have not been provided by the model. We use here
875  * a simple method as this method is likely used for spatial model parameters,
876  * and the spatial model may eventually be noisy due to numerical integration
877  * limits.
878  ***************************************************************************/
879 double GObservation::model_grad(const GModel& model,
880  const GModelPar& par,
881  const GEvent& event) const
882 {
883  // Initialise gradient
884  double grad = 0.0;
885 
886  // Compute gradient only if parameter is free
887  if (par.is_free()) {
888 
889  // Get non-const model pointer
890  GModelPar* ptr = const_cast<GModelPar*>(&par);
891 
892  // Save current model parameter
893  GModelPar current = par;
894 
895  // Get actual parameter value
896  double x = par.factor_value();
897 
898  // Set step size for computation of derivative.
899  double h = m_grad_step_size;
900 
901  // Re-adjust the step-size h in case that the initial step size is
902  // larger than the allowed parameter range
903  if (par.has_min() && par.has_max()) {
904  double par_h = par.factor_max() - par.factor_min();
905  if (par_h < h) {
906  h = par_h;
907  }
908  }
909 
910  // Continue only if step size is positive
911  if (h > 0.0) {
912 
913  // Remove any boundaries to avoid limitations
914  //ptr->remove_range(); // Not needed in principle
915 
916  // Setup derivative function
917  GObservation::model_func function(this, &model, ptr, &event);
918  GDerivative derivative(&function);
919 
920  // If we are too close to the minimum boundary use a right sided
921  // difference ...
922  if (par.has_min() && ((x-par.factor_min()) < h)) {
923  grad = derivative.right_difference(x, h);
924  }
925 
926  // ... otherwise if we are too close to the maximum boundary use
927  // a left sided difference ...
928  else if (par.has_max() && ((par.factor_max()-x) < h)) {
929  grad = derivative.left_difference(x, h);
930  }
931 
932  // ... otherwise use a symmetric difference
933  else {
934  grad = derivative.difference(x, h);
935  }
936 
937  } // endif: step size was positive
938 
939  // Restore current model parameter
940  *ptr = current;
941 
942  } // endif: model parameter was free
943 
944  // Return gradient
945  return grad;
946 }
947 
948 
949 /***********************************************************************//**
950  * @brief Returns parameter gradients of model for all events
951  *
952  * @param[in] model Model.
953  * @param[in] par Model parameter.
954  *
955  * This method uses a robust but simple difference method to estimate
956  * parameter gradients that have not been provided by the model. We use here
957  * a simple method as this method is likely used for spatial model parameters,
958  * and the spatial model may eventually be noisy due to numerical integration
959  * limits.
960  ***************************************************************************/
962  const GModelPar& par) const
963 {
964  // Initialise gradients
965  GVector gradients(events()->size());
966 
967  // Compute gradient only if parameter is free
968  if (par.is_free()) {
969 
970  // Get non-const model pointer
971  GModelPar* ptr = const_cast<GModelPar*>(&par);
972 
973  // Save current model parameter
974  GModelPar current = par;
975 
976  // Get actual parameter value
977  double x = par.factor_value();
978 
979  // Set step size for computation of derivative.
980  double h = m_grad_step_size;
981 
982  // Re-adjust the step-size h in case that the initial step size is
983  // larger than the allowed parameter range
984  if (par.has_min() && par.has_max()) {
985  double par_h = par.factor_max() - par.factor_min();
986  if (par_h < h) {
987  h = par_h;
988  }
989  }
990 
991  // Continue only if step size is positive
992  if (h > 0.0) {
993 
994  // Setup derivative function
995  //GObservation::model_func function(this, &model, ptr, &event);
996  //GDerivative derivative(&function);
997 
998  // If we are too close to the minimum boundary use a right sided
999  // difference ...
1000  if (par.has_min() && ((x-par.factor_min()) < h)) {
1001  //grad = derivative.right_difference(x, h);
1002  // Compute fs1
1003  ptr->factor_value(x);
1004  GVector fs1 = model.eval(*this, NULL);
1005 
1006  // Compute fs2
1007  ptr->factor_value(x+h);
1008  GVector fs2 = model.eval(*this, NULL);
1009 
1010  // Compute derivative
1011  gradients = (fs1 - fs2) / h;
1012 
1013  }
1014 
1015  // ... otherwise if we are too close to the maximum boundary use
1016  // a left sided difference ...
1017  else if (par.has_max() && ((par.factor_max()-x) < h)) {
1018  //grad = derivative.left_difference(x, h);
1019  // Compute fs1
1020  ptr->factor_value(x);
1021  GVector fs1 = model.eval(*this, NULL);
1022 
1023  // Compute fs2
1024  ptr->factor_value(x-h);
1025  GVector fs2 = model.eval(*this, NULL);
1026 
1027  // Compute derivative
1028  gradients = (fs1 - fs2) / h;
1029  }
1030 
1031  // ... otherwise use a symmetric difference
1032  else {
1033  //grad = derivative.difference(x, h);
1034  // Compute fs1
1035  ptr->factor_value(x+h);
1036  GVector fs1 = model.eval(*this, NULL);
1037 
1038  // Compute fs2
1039  ptr->factor_value(x-h);
1040  GVector fs2 = model.eval(*this, NULL);
1041 
1042  // Compute derivative
1043  gradients = 0.5 * (fs1 - fs2) / h;
1044  }
1045 
1046  } // endif: step size was positive
1047 
1048  // Restore current model parameter
1049  *ptr = current;
1050 
1051  } // endif: model parameter was free
1052 
1053  // Return gradients
1054  return gradients;
1055 }
1056 
1057 
1058 /***********************************************************************//**
1059  * @brief Returns parameter gradient of Npred
1060  *
1061  * @param[in] model Gamma-ray source model.
1062  * @param[in] par Model parameter.
1063  *
1064  * Computes
1065  *
1066  * \f[
1067  * \frac{{\rm d} N_{\rm pred}}{{\rm d} a_i}
1068  * \f]
1069  *
1070  * where
1071  *
1072  * \f[
1073  * N_{\rm pred} = \int_{\rm GTI} \int_{E_{\rm bounds}} \int_{\rm ROI}
1074  * P(p',E',t') \, dp' \, dE' \, dt'
1075  * \f]
1076  *
1077  * is the integral of the event probability
1078  *
1079  * \f[
1080  * P(p',E',t') = \int \int \int
1081  * S(p,E,t) \times R(p',E',t'|p,E,t) \, dp \, dE \, dt
1082  * \f]
1083  *
1084  * and
1085  * \f$a_i\f$ is the model parameter \f$i\f$.
1086  * Furthermore
1087  * \f$S(p,E,t)\f$ is the source model,
1088  * \f$R(p',E',t'|p,E,t)\f$ is the instrument response function,
1089  * \f$p'\f$ is the measured photon direction,
1090  * \f$E'\f$ is the measured photon energy,
1091  * \f$t'\f$ is the measured photon arrival time,
1092  * \f$p\f$ is the true photon arrival direction,
1093  * \f$E\f$ is the true photon energy, and
1094  * \f$t\f$ is the true photon arrival time.
1095  *
1096  * This method uses a robust but simple difference method to estimate
1097  * parameter gradients that have not been provided by the model. We use here
1098  * a simple method as this method is likely used for spatial model parameters,
1099  * and the spatial model may eventually be noisy due to numerical integration
1100  * limits.
1101  ***************************************************************************/
1102 double GObservation::npred_grad(const GModel& model, const GModelPar& par) const
1103 {
1104  // Initialise result
1105  double grad = 0.0;
1106 
1107  // Compute gradient only if parameter is free
1108  if (par.is_free()) {
1109 
1110  // Get non-const model pointer
1111  GModelPar* ptr = const_cast<GModelPar*>(&par);
1112 
1113  // Save current model parameter
1114  GModelPar current = par;
1115 
1116  // Get actual parameter value
1117  double x = par.factor_value();
1118 
1119  // Set step size for computation of derivative.
1120  double h = m_grad_step_size;
1121 
1122  // Re-adjust the step-size h in case that the initial step size is
1123  // larger than the allowed parameter range
1124  if (par.has_min() && par.has_max()) {
1125  double par_h = par.factor_max() - par.factor_min();
1126  if (par_h < h) {
1127  h = par_h;
1128  }
1129  }
1130 
1131  // Continue only if step size is positive
1132  if (h > 0.0) {
1133 
1134  // Remove any boundaries to avoid limitations
1135  //ptr->remove_range(); // Not needed in principle
1136 
1137  // Setup derivative function
1138  GObservation::npred_func function(this, &model, ptr);
1139  GDerivative derivative(&function);
1140 
1141  // If we are too close to the minimum boundary use a right sided
1142  // difference ...
1143  if (par.has_min() && ((x-par.factor_min()) < h)) {
1144  grad = derivative.right_difference(x, h);
1145  }
1146 
1147  // ... otherwise if we are too close to the maximum boundary use
1148  // a left sided difference ...
1149  else if (par.has_max() && ((par.factor_max()-x) < h)) {
1150  grad = derivative.left_difference(x, h);
1151  }
1152 
1153  // ... otherwise use a symmetric difference
1154  else {
1155  grad = derivative.difference(x, h);
1156  }
1157 
1158  } // endif: step size was positive
1159 
1160  // Restore current model parameter
1161  *ptr = current;
1162 
1163  } // endif: model parameter was free
1164 
1165  // Return result
1166  return grad;
1167 }
1168 
1169 
1170 /***********************************************************************//**
1171  * @brief Response cache removal hook
1172  *
1173  * @param[in] name Model name.
1174  *
1175  * Remove response cache for model @p name from response cache.
1176  ***************************************************************************/
1177 void GObservation::remove_response_cache(const std::string& name)
1178 {
1179  // Build model name
1180  std::string model_name = id() + ":" + name;
1181 
1182  // Remove response cache
1183  const_cast<GResponse*>(this->response())->remove_response_cache(model_name);
1184 
1185  // Return
1186  return;
1187 }
1188 
1189 
1190 /***********************************************************************//**
1191  * @brief Check whether a model parameter has an analytical gradient
1192  *
1193  * @param[in] model Model.
1194  * @param[in] par Model parameter.
1195  * @returns True if model parameter is free and has an analytical gradient.
1196  *
1197  * Checks whether a model parameter is free and has an analytical gradient.
1198  ***************************************************************************/
1199 bool GObservation::has_gradient(const GModel& model, const GModelPar& par) const
1200 {
1201  // Initialise flag
1202  bool has_gradient = false;
1203 
1204  // Only consider free parameters
1205  if (par.is_free()) {
1206 
1207  // Build identifier
1208  std::string id = model.name() + ":" + par.name();
1209 
1210  // Search for parameter address in stack and set flag to true if
1211  // parameter address was found
1212  for (int i = 0; i < m_pars_with_gradients.size(); ++i) {
1213  if (m_pars_with_gradients[i] == id) {
1214  has_gradient = true;
1215  break;
1216  }
1217  }
1218 
1219  } // endif: parameter was free
1220 
1221  // Return flag
1222  return has_gradient;
1223 }
1224 
1225 
1226 /***********************************************************************//**
1227  * @brief Signals that an analytical gradient was computed for a model
1228  * parameter
1229  *
1230  * @param[in] model Model.
1231  * @param[in] par Model parameter.
1232  *
1233  * Signals that an analytical gradient was computed for a model parameter.
1234  ***************************************************************************/
1235 void GObservation::computed_gradient(const GModel& model, const GModelPar& par) const
1236 {
1237  // Initialise flag
1238  bool in_stack = false;
1239 
1240  // Build identifier
1241  std::string id = model.name() + ":" + par.name();
1242 
1243  // Check if parameter is already in stack
1244  for (int i = 0; i < m_pars_with_gradients.size(); ++i) {
1245  if (m_pars_with_gradients[i] == id) {
1246  in_stack = true;
1247  break;
1248  }
1249  }
1250 
1251  // If parameter is not yet in stack the push it on the stack
1252  if (!in_stack) {
1253  m_pars_with_gradients.push_back(id);
1254  }
1255 
1256  // Return
1257  return;
1258 }
1259 
1260 
1261 /*==========================================================================
1262  = =
1263  = Private methods =
1264  = =
1265  ==========================================================================*/
1266 
1267 /***********************************************************************//**
1268  * @brief Initialise class members
1269  *
1270  * The step size for gradient computation is fixed by default to 0.0002
1271  * degrees, which corresponds to about 1 arcsec for parameters that are
1272  * given in degrees. The reasoning behind this value is that parameters that
1273  * use numerical gradients are typically angles, such as for example the
1274  * position, and we want to achieve arcsec precision with this method.
1275  ***************************************************************************/
1277 {
1278  // Initialise members
1279  m_name.clear();
1280  m_id.clear();
1281  m_statistic = "cstat";
1282  m_events = NULL;
1283  m_grad_step_size = 0.0002;
1284 
1285  // Initialise stack of identifiers of parameters with gradients
1286  m_pars_with_gradients.clear();
1287 
1288  // Return
1289  return;
1290 }
1291 
1292 
1293 /***********************************************************************//**
1294  * @brief Copy class members
1295  *
1296  * @param[in] obs Observation.
1297  *
1298  * Copy members from an observation.
1299  ***************************************************************************/
1301 {
1302  // Copy members
1303  m_name = obs.m_name;
1304  m_id = obs.m_id;
1305  m_statistic = obs.m_statistic;
1307 
1308  // Copy stack of identifiers of parameters with gradients
1310 
1311  // Clone members
1312  m_events = (obs.m_events != NULL) ? obs.m_events->clone() : NULL;
1313 
1314  // Return
1315  return;
1316 }
1317 
1318 
1319 /***********************************************************************//**
1320  * @brief Delete class members
1321  ***************************************************************************/
1323 {
1324  // Free members
1325  if (m_events != NULL) delete m_events;
1326 
1327  // Signal free pointers
1328  m_events = NULL;
1329 
1330  // Return
1331  return;
1332 }
1333 
1334 
1335 /*==========================================================================
1336  = =
1337  = Likelihood methods =
1338  = =
1339  ==========================================================================*/
1340 
1341 /***********************************************************************//**
1342  * @brief Evaluate log-likelihood function for Poisson statistic and
1343  * unbinned analysis (version with working arrays)
1344  *
1345  * @param[in] models Models.
1346  * @param[in,out] gradients Gradient.
1347  * @param[in,out] curvature Curvature matrix.
1348  * @param[in,out] npred Number of predicted events.
1349  * @return Likelihood value.
1350  *
1351  * This method evaluates the -(log-likelihood) function for parameter
1352  * optimisation using unbinned analysis and Poisson statistic.
1353  * The -(log-likelihood) function is given by
1354  *
1355  * \f[
1356  * L = N_{\rm pred} - \sum_i \log e_i
1357  * \f]
1358  *
1359  * where
1360  * \f$N_{\rm pred}\f$ is the number of events predicted by the model, and
1361  * the sum is taken over all events. This method also computes the
1362  * parameter gradients
1363  * \f$\delta L/dp\f$
1364  * and the curvature matrix
1365  * \f$\delta^2 L/dp_1 dp_2\f$.
1366  ***************************************************************************/
1368  GVector* gradients,
1369  GMatrixSparse* curvature,
1370  double* npred) const
1371 {
1372  // Initialise likelihood value
1373  double value = 0.0;
1374 
1375  // Get number of events and parameters
1376  int nevents = events()->size();
1377  int npars = gradients->size();
1378 
1379  // Allocate some working arrays
1380  int* inx = new int[npars];
1381  double* values = new double[npars];
1382  GMatrixSparse wrk_matrix(nevents, npars);
1383  GVector wrk_grad(npars);
1384 
1385  // Determine Npred value and gradient for this observation
1386  double npred_value = this->npred(models, &wrk_grad);
1387 
1388  // Update likelihood, Npred and gradients
1389  value += npred_value;
1390  *npred += npred_value;
1391  *gradients += wrk_grad;
1392 
1393  // Compute model and derivative
1394  GVector model_vector = this->model(models, &wrk_matrix);
1395 
1396  // Transpose matrix
1397  wrk_matrix = wrk_matrix.transpose();
1398 
1399  // Iterate over all events
1400  for (int i = 0; i < nevents; ++i) {
1401 
1402  // Skip events that should not be used
1403  if (!use_event_for_likelihood(i)) {
1404  continue;
1405  }
1406 
1407  // Get model value
1408  double model = model_vector[i];
1409 
1410  // Skip bin if model is too small (avoids -Inf or NaN gradients)
1411  if (model <= minmod) {
1412  continue;
1413  }
1414 
1415  // Extract working gradient multiplied by bin size
1416  GVector wrk_grad = wrk_matrix.column(i);
1417 
1418  // Create index array of non-zero derivatives and initialise working
1419  // array
1420  int ndev = 0;
1421  for (int ipar = 0; ipar < npars; ++ipar) {
1422  values[ipar] = 0.0;
1423  if (wrk_grad[ipar] != 0.0 && !gammalib::is_infinite(wrk_grad[ipar])) {
1424  inx[ndev] = ipar;
1425  ndev++;
1426  }
1427  }
1428 
1429  // Update Poissonian statistic (excluding factorial term for faster
1430  // computation)
1431  value -= std::log(model);
1432 
1433  // Skip bin now if there are no non-zero derivatives
1434  if (ndev < 1) {
1435  continue;
1436  }
1437 
1438  // Update gradient vector and curvature matrix.
1439  double fb = 1.0 / model;
1440  double fa = fb / model;
1441  for (int jdev = 0; jdev < ndev; ++jdev) {
1442 
1443  // Initialise computation
1444  register int jpar = inx[jdev];
1445  double g = wrk_grad[jpar];
1446  double fa_i = fa * g;
1447 
1448  // Update gradient.
1449  (*gradients)[jpar] -= fb * g;
1450 
1451  // Loop over rows
1452  register int* ipar = inx;
1453 
1454  for (register int idev = 0; idev < ndev; ++idev, ++ipar) {
1455  values[idev] = fa_i * wrk_grad[*ipar];
1456  }
1457 
1458  // Add column to matrix
1459  curvature->add_to_column(jpar, values, inx, ndev);
1460 
1461  } // endfor: looped over columns
1462 
1463  } // endfor: iterated over all events
1464 
1465  // Free temporary memory
1466  if (values != NULL) delete [] values;
1467  if (inx != NULL) delete [] inx;
1468 
1469  // Return
1470  return value;
1471 }
1472 
1473 
1474 /***********************************************************************//**
1475  * @brief Evaluate log-likelihood function for Poisson statistic and
1476  * binned analysis (version with working arrays)
1477  *
1478  * @param[in] models Models.
1479  * @param[in,out] gradients Gradient.
1480  * @param[in,out] curvature Curvature matrix.
1481  * @param[in,out] npred Number of predicted events.
1482  * @return Likelihood value.
1483  *
1484  * This method evaluates the -(log-likelihood) function for parameter
1485  * optimisation using binned analysis and Poisson statistic.
1486  * The -(log-likelihood) function is given by
1487  *
1488  * \f[
1489  * L=-\sum_i n_i \log e_i - e_i
1490  * \f]
1491  *
1492  * where the sum is taken over all data space bins, \f$n_i\f$ is the
1493  * observed number of counts and \f$e_i\f$ is the model.
1494  * This method also computes the parameter gradients
1495  * \f$\delta L/dp\f$
1496  * and the curvature matrix
1497  * \f$\delta^2 L/dp_1 dp_2\f$
1498  * and also updates the total number of predicted events m_npred.
1499  ***************************************************************************/
1501  GVector* gradients,
1502  GMatrixSparse* curvature,
1503  double* npred) const
1504 {
1505  // Initialise likelihood value
1506  double value = 0.0;
1507 
1508  // Initialise statistic
1509  #if defined(G_OPT_DEBUG)
1510  int n_bins = 0;
1511  int n_used = 0;
1512  int n_small_model = 0;
1513  int n_zero_data = 0;
1514  int n_exclude_data = 0;
1515  double sum_data = 0.0;
1516  double sum_model = 0.0;
1517  double init_value = value;
1518  #endif
1519 
1520  // Get number of events and parameters
1521  int nevents = events()->size();
1522  int npars = gradients->size();
1523 
1524  // Allocate some working arrays
1525  int* inx = new int[npars];
1526  double* values = new double[npars];
1527  GMatrixSparse wrk_matrix(nevents, npars);
1528 
1529  // Compute model and derivative
1530  GVector model_vector = this->model(models, &wrk_matrix);
1531 
1532  // Transpose matrix
1533  wrk_matrix = wrk_matrix.transpose();
1534 
1535  // Iterate over all bins
1536  for (int i = 0; i < nevents; ++i) {
1537 
1538  // Skip events that should not be used
1539  if (!use_event_for_likelihood(i)) {
1540  continue;
1541  }
1542 
1543  // Update number of bins
1544  #if defined(G_OPT_DEBUG)
1545  n_bins++;
1546  #endif
1547 
1548  // Get event pointer
1549  const GEventBin* bin =
1550  (*(static_cast<GEventCube*>(const_cast<GEvents*>(events()))))[i];
1551 
1552  // Get number of counts in bin
1553  double data = bin->counts();
1554 
1555  // Skip bin if data is negative (filtering flag)
1556  if (data < 0) {
1557  #if defined(G_OPT_DEBUG)
1558  n_exclude_data++;
1559  #endif
1560  continue;
1561  }
1562 
1563  // Get model value
1564  double model = model_vector[i] * bin->size();
1565 
1566  // Skip bin if model is too small (avoids -Inf or NaN gradients)
1567  if (model <= minmod) {
1568  #if defined(G_OPT_DEBUG)
1569  n_small_model++;
1570  #endif
1571  continue;
1572  }
1573 
1574  // Update statistic
1575  #if defined(G_OPT_DEBUG)
1576  n_used++;
1577  sum_data += data;
1578  sum_model += model;
1579  #endif
1580 
1581  // Update Npred
1582  *npred += model;
1583 
1584  // Extract working gradient multiplied by bin size
1585  GVector wrk_grad = wrk_matrix.column(i) * bin->size();
1586 
1587  // Create index array of non-zero derivatives and initialise working
1588  // array
1589  int ndev = 0;
1590  for (int ipar = 0; ipar < npars; ++ipar) {
1591  values[ipar] = 0.0;
1592  if (wrk_grad[ipar] != 0.0 && !gammalib::is_infinite(wrk_grad[ipar])) {
1593  inx[ndev] = ipar;
1594  ndev++;
1595  }
1596  }
1597 
1598  // Update gradient vector and curvature matrix. To avoid
1599  // unneccessary computations we distinguish the case where
1600  // data>0 and data=0. The second case requires much less
1601  // computation since it does not contribute to the curvature
1602  // matrix ...
1603  if (data > 0.0) {
1604 
1605  // Update Poissonian statistic (excluding factorial term for
1606  // faster computation)
1607  value -= data * std::log(model) - model;
1608 
1609  // Skip bin now if there are no non-zero derivatives
1610  if (ndev < 1) {
1611  continue;
1612  }
1613 
1614  // Pre computation
1615  double fb = data / model;
1616  double fc = (1.0 - fb);
1617  double fa = fb / model;
1618 
1619  // Loop over columns
1620  for (int jdev = 0; jdev < ndev; ++jdev) {
1621 
1622  // Initialise computation
1623  register int jpar = inx[jdev];
1624  double g = wrk_grad[jpar];
1625  double fa_i = fa * g;
1626 
1627  // Update gradient
1628  (*gradients)[jpar] += fc * g;
1629 
1630  // Loop over rows
1631  register int* ipar = inx;
1632  for (register int idev = 0; idev < ndev; ++idev, ++ipar) {
1633  values[idev] = fa_i * wrk_grad[*ipar];
1634  }
1635 
1636  // Add column to matrix
1637  curvature->add_to_column(jpar, values, inx, ndev);
1638 
1639  } // endfor: looped over columns
1640 
1641  } // endif: data was > 0
1642 
1643  // ... handle now data=0
1644  else {
1645 
1646  // Update statistic
1647  #if defined(G_OPT_DEBUG)
1648  n_zero_data++;
1649  #endif
1650 
1651  // Update Poissonian statistic (excluding factorial term for
1652  // faster computation)
1653  value += model;
1654 
1655  // Skip bin now if there are no non-zero derivatives
1656  if (ndev < 1) {
1657  continue;
1658  }
1659 
1660  // Update gradient
1661  register int* ipar = inx;
1662  for (register int idev = 0; idev < ndev; ++idev, ++ipar) {
1663  (*gradients)[*ipar] += wrk_grad[*ipar];
1664  }
1665 
1666  } // endif: data was 0
1667 
1668  } // endfor: iterated over all events
1669 
1670  // Free temporary memory
1671  if (values != NULL) delete [] values;
1672  if (inx != NULL) delete [] inx;
1673 
1674  // Dump statistic
1675  #if defined(G_OPT_DEBUG)
1676  std::cout << "Number of bins: " << n_bins << std::endl;
1677  std::cout << "Number of bins used for computation: " << n_used << std::endl;
1678  std::cout << "Number of bins excluded from data: " << n_exclude_data << std::endl;
1679  std::cout << "Number of bins excluded due to small model: " << n_small_model << std::endl;
1680  std::cout << "Number of bins with zero data: " << n_zero_data << std::endl;
1681  std::cout << "Sum of data: " << sum_data << std::endl;
1682  std::cout << "Sum of model: " << sum_model << std::endl;
1683  std::cout << "Initial statistic: " << init_value << std::endl;
1684  std::cout << "Statistic: " << value-init_value << std::endl;
1685  #endif
1686 
1687  // Return
1688  return value;
1689 }
1690 
1691 
1692 /***********************************************************************//**
1693  * @brief Evaluate log-likelihood function for Gaussian statistic and
1694  * binned analysis (version with working arrays)
1695  *
1696  * @param[in] models Models.
1697  * @param[in,out] gradients Gradient.
1698  * @param[in,out] curvature Curvature matrix.
1699  * @param[in,out] npred Number of predicted events.
1700  * @return Likelihood value.
1701  *
1702  * This method evaluates the -(log-likelihood) function for parameter
1703  * optimisation using binned analysis and Poisson statistic.
1704  * The -(log-likelihood) function is given by
1705  *
1706  * \f[
1707  * L = 1/2 \sum_i (n_i - e_i)^2 \sigma_i^{-2}
1708  * \f]
1709  *
1710  * where the sum is taken over all data space bins, \f$n_i\f$ is the
1711  * observed number of counts, \f$e_i\f$ is the model and \f$\sigma_i\f$
1712  * is the statistical uncertainty.
1713  * This method also computes the parameter gradients
1714  * \f$\delta L/dp\f$
1715  * and the curvature matrix
1716  * \f$\delta^2 L/dp_1 dp_2\f$
1717  * and also updates the total number of predicted events m_npred.
1718  ***************************************************************************/
1720  GVector* gradients,
1721  GMatrixSparse* curvature,
1722  double* npred) const
1723 {
1724  // Initialise likelihood value
1725  double value = 0.0;
1726 
1727  // Get number of events and parameters
1728  int nevents = events()->size();
1729  int npars = gradients->size();
1730 
1731  // Allocate some working arrays
1732  int* inx = new int[npars];
1733  double* values = new double[npars];
1734  GMatrixSparse wrk_matrix(nevents, npars);
1735 
1736  // Compute model and derivative
1737  GVector model_vector = this->model(models, &wrk_matrix);
1738 
1739  // Transpose matrix
1740  wrk_matrix = wrk_matrix.transpose();
1741 
1742  // Iterate over all bins
1743  for (int i = 0; i < nevents; ++i) {
1744 
1745  // Skip events that should not be used
1746  if (!use_event_for_likelihood(i)) {
1747  continue;
1748  }
1749 
1750  // Get event pointer
1751  const GEventBin* bin =
1752  (*(static_cast<GEventCube*>(const_cast<GEvents*>(events()))))[i];
1753 
1754  // Get number of counts in bin
1755  double data = bin->counts();
1756 
1757  // Skip bin if data is negative (filtering flag)
1758  if (data < 0) {
1759  continue;
1760  }
1761 
1762  // Get statistical uncertainty
1763  double sigma = bin->error();
1764 
1765  // Skip bin if statistical uncertainty is too small
1766  if (sigma <= minerr) {
1767  continue;
1768  }
1769 
1770  // Get model value
1771  double model = model_vector[i] * bin->size();
1772 
1773  // Skip bin if model is too small (avoids -Inf or NaN gradients)
1774  if (model <= minmod) {
1775  continue;
1776  }
1777 
1778  // Update Npred
1779  *npred += model;
1780 
1781  // Extract working gradient multiplied by bin size
1782  GVector wrk_grad = wrk_matrix.column(i) * bin->size();
1783 
1784  // Create index array of non-zero derivatives and initialise working
1785  // array
1786  int ndev = 0;
1787  for (int ipar = 0; ipar < npars; ++ipar) {
1788  values[ipar] = 0.0;
1789  if (wrk_grad[ipar] != 0.0 && !gammalib::is_infinite(wrk_grad[ipar])) {
1790  inx[ndev] = ipar;
1791  ndev++;
1792  }
1793  }
1794 
1795  // Set weight
1796  double weight = 1.0 / (sigma * sigma);
1797 
1798  // Update Gaussian statistic
1799  double fa = data - model;
1800  value += 0.5 * (fa * fa * weight);
1801 
1802  // Skip bin now if there are no non-zero derivatives
1803  if (ndev < 1) {
1804  continue;
1805  }
1806 
1807  // Loop over columns
1808  for (int jdev = 0; jdev < ndev; ++jdev) {
1809 
1810  // Initialise computation
1811  register int jpar = inx[jdev];
1812  double fa_i = wrk_grad[jpar] * weight;
1813 
1814  // Update gradient
1815  (*gradients)[jpar] -= fa * fa_i;
1816 
1817  // Loop over rows
1818  register int* ipar = inx;
1819  for (register int idev = 0; idev < ndev; ++idev, ++ipar) {
1820  values[idev] = fa_i * wrk_grad[*ipar];
1821  }
1822 
1823  // Add column to matrix
1824  curvature->add_to_column(jpar, values, inx, ndev);
1825 
1826  } // endfor: looped over columns
1827 
1828  } // endfor: iterated over all events
1829 
1830  // Free temporary memory
1831  if (values != NULL) delete [] values;
1832  if (inx != NULL) delete [] inx;
1833 
1834  // Return
1835  return value;
1836 }
1837 
1838 
1839 /***********************************************************************//**
1840  * @brief Check whether bin should be used for likelihood analysis
1841  *
1842  * @param[in] index Event index.
1843  * @return True.
1844  *
1845  * This is a dummy virtual method that allows implementation of a hook for
1846  * event selection in the likelihood computation. The dummy method always
1847  * returns true.
1848  ***************************************************************************/
1849 bool GObservation::use_event_for_likelihood(const int& index) const
1850 {
1851  // Return true
1852  return true;
1853 }
1854 
1855 
1856 /*==========================================================================
1857  = =
1858  = Model gradient methods =
1859  = =
1860  ==========================================================================*/
1861 
1862 /***********************************************************************//**
1863  * @brief Model function evaluation for gradient computation
1864  *
1865  * @param[in] x Function value.
1866  ***************************************************************************/
1867 double GObservation::model_func::eval(const double& x)
1868 {
1869  // Set value
1870  m_par->factor_value(x);
1871 
1872  // Compute model value
1873  double value = m_model->eval(*m_event, *m_parent);
1874 
1875  // Return value
1876  return value;
1877 }
1878 
1879 
1880 /*==========================================================================
1881  = =
1882  = Npred computation methods =
1883  = =
1884  ==========================================================================*/
1885 
1886 /***********************************************************************//**
1887  * @brief Npred function evaluation for gradient computation
1888  *
1889  * @param[in] x Function value.
1890  ***************************************************************************/
1891 double GObservation::npred_func::eval(const double& x)
1892 {
1893  // Set value
1894  m_par->factor_value(x);
1895 
1896  // Compute Npred value
1897  double npred = m_parent->npred(*m_model);
1898 
1899  // Return value
1900  return npred;
1901 }
1902 
1903 
1904 /***********************************************************************//**
1905  * @brief Integration kernel for npred() method
1906  *
1907  * @param[in] x Function value.
1908  ***************************************************************************/
1909 double GObservation::npred_kern::eval(const double& x)
1910 {
1911  // Convert argument in native reference in seconds
1912  GTime time;
1913  time.secs(x);
1914 
1915  // Return value
1916  return (m_parent->npred_spec(*m_model, time));
1917 }
1918 
1919 
1920 /***********************************************************************//**
1921  * @brief Integrates spatially integrated Npred kernel spectrally
1922  *
1923  * @param[in] model Gamma-ray source model.
1924  * @param[in] obsTime Measured photon arrival time.
1925  *
1926  * @exception GException::erange_invalid
1927  * Energy range is invalid.
1928  *
1929  * Computes
1930  *
1931  * \f[
1932  * N_{\rm pred}(t') = \int_{E_{\rm bounds}} \int_{\rm ROI}
1933  * P(p',E',t') \, dp' \, dE'
1934  * \f]
1935  *
1936  * of the event probability
1937  *
1938  * \f[
1939  * P(p',E',t') = \int \int \int
1940  * S(p,E,t) \times R(p',E',t'|p,E,t) \, dp \, dE \, dt
1941  * \f]
1942  *
1943  * where
1944  * \f$S(p,E,t)\f$ is the source model,
1945  * \f$R(p',E',t'|p,E,t)\f$ is the instrument response function,
1946  * \f$p'\f$ is the measured photon direction,
1947  * \f$E'\f$ is the measured photon energy,
1948  * \f$t'\f$ is the measured photon arrival time,
1949  * \f$p\f$ is the true photon arrival direction,
1950  * \f$E\f$ is the true photon energy, and
1951  * \f$t\f$ is the true photon arrival time.
1952  *
1953  * This method performs the energy intergration over the energy boundaries
1954  * \f$E_{\rm bounds}\f$.
1955  ***************************************************************************/
1957  const GTime& obsTime) const
1958 {
1959  // Set number of iterations for Romberg integration.
1960  static const int iter = 8;
1961 
1962  // Initialise result
1963  double result = 0.0;
1964 
1965  // Get energy boundaries
1966  GEbounds ebounds = events()->ebounds();
1967 
1968  // Loop over energy boundaries
1969  for (int i = 0; i < ebounds.size(); ++i) {
1970 
1971  // Get boundaries in MeV
1972  double emin = ebounds.emin(i).MeV();
1973  double emax = ebounds.emax(i).MeV();
1974 
1975  // Continue only if valid
1976  if (emax > emin) {
1977 
1978  // Setup integration function
1979  GObservation::npred_spec_kern integrand(this, &model, &obsTime);
1980  GIntegral integral(&integrand);
1981 
1982  // Set number of iterations
1983  integral.fixed_iter(iter);
1984 
1985  // Do Romberg integration
1986  #if defined(G_LN_ENERGY_INT)
1987  emin = std::log(emin);
1988  emax = std::log(emax);
1989  #endif
1990  result += integral.romberg(emin, emax);
1991 
1992  } // endif: energy interval was valid
1993 
1994  } // endfor: looped over energy boundaries
1995 
1996  // Compile option: Check for NaN/Inf
1997  #if defined(G_NAN_CHECK)
1998  if (gammalib::is_notanumber(result) || gammalib::is_infinite(result)) {
1999  std::cout << "*** ERROR: GObservation::npred_spec:";
2000  std::cout << " NaN/Inf encountered";
2001  std::cout << " (result=" << result;
2002  std::cout << ")" << std::endl;
2003  }
2004  #endif
2005 
2006  // Return result
2007  return result;
2008 }
2009 
2010 
2011 /***********************************************************************//**
2012  * @brief Integration kernel for npred_spec() method
2013  *
2014  * @param[in] x Function value.
2015  *
2016  * This method implements the integration kernel needed for the npred_spec()
2017  * method. If G_LN_ENERGY_INT is defined the energy integration is done
2018  * logarithmically, i.e. @p x is given in ln(energy) instead of energy.
2019  ***************************************************************************/
2021 {
2022  // Set energy
2023  GEnergy eng;
2024  #if defined(G_LN_ENERGY_INT)
2025  double expx = std::exp(x);
2026  eng.MeV(expx);
2027  #else
2028  eng.MeV(x);
2029  #endif
2030 
2031  // Get function value
2032  double value = m_model->npred(eng, *m_time, *m_parent);
2033 
2034  // Save value if needed
2035  #if defined(G_NAN_CHECK)
2036  double value_out = value;
2037  #endif
2038 
2039  // Correct for variable substitution
2040  #if defined(G_LN_ENERGY_INT)
2041  value *= expx;
2042  #endif
2043 
2044  // Compile option: Check for NaN
2045  #if defined(G_NAN_CHECK)
2046  if (gammalib::is_notanumber(value) || gammalib::is_infinite(value)) {
2047  std::cout << "*** ERROR: GObservation::npred_spec_kern::eval";
2048  std::cout << "(x=" << x << "): ";
2049  std::cout << " NaN/Inf encountered";
2050  std::cout << " (value=" << value;
2051  std::cout << " (value_out=" << value_out;
2052  #if defined(G_LN_ENERGY_INT)
2053  std::cout << " exp(x)=" << expx;
2054  #endif
2055  std::cout << ")" << std::endl;
2056  }
2057  #endif
2058 
2059  // Return value
2060  return value;
2061 }
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:100
double eval(const double &x)
Integration kernel for npred_spec() method.
virtual bool is_constant(void) const =0
virtual double likelihood_gaussian_binned(const GModels &models, GVector *gradients, GMatrixSparse *curvature, double *npred) const
Evaluate log-likelihood function for Gaussian statistic and binned analysis (version with working arr...
double m_grad_step_size
Gradient step size.
const GEvent * m_event
Event.
void computed_gradient(const GModel &model, const GModelPar &par) const
Signals that an analytical gradient was computed for a model parameter.
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 factor gradient.
virtual double ontime(void) const =0
double romberg(std::vector< double > bounds, const int &order=5)
Perform Romberg integration.
Definition: GIntegral.cpp:381
std::string m_name
Observation name.
const std::string & name(void) const
Return parameter name.
Sparse matrix class interface definition.
virtual std::string instrument(void) const =0
Model container class definition.
#define G_MODEL2
void ebounds(const GEbounds &ebounds)
Set energy boundaries.
Definition: GEvents.cpp:136
int size(void) const
Return number of energy boundaries.
Definition: GEbounds.hpp:163
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
virtual double npred(const GModels &models, GVector *gradients=NULL) const
Return total number (and optionally gradients) of predicted counts for all models.
virtual double likelihood(const GModels &models, GVector *gradients, GMatrixSparse *curvature, double *npred) const
Compute likelihood function.
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
const bool & has_eval_indices(void) const
Signals that parameter indices updated by eval() method were set.
Definition: GModel.hpp:408
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:55
Gammalib tools definition.
GIntegral class interface definition.
Definition: GIntegral.hpp:46
const std::vector< int > & eval_indices(void) const
Return vector of parameter indices updated by eval() method.
Definition: GModel.hpp:420
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.
GMatrixSparse transpose(void) const
Return transposed matrix.
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:201
Model parameter class interface definition.
int size(void) const
Return number of parameters in model.
Definition: GModel.hpp:233
bool is_infinite(const double &x)
Signal if argument is infinite.
Definition: GTools.hpp:184
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:570
#define G_MODEL1
Model parameter class.
Definition: GModelPar.hpp:87
Model container class.
Definition: GModels.hpp:152
#define G_NPRED
const std::string & id(void) const
Return observation identifier.
Energy boundaries container class.
Definition: GEbounds.hpp:60
const std::string & name(void) const
Return parameter name.
Definition: GModel.hpp:261
const double & factor_max(void) const
Return parameter maximum factor boundary.
virtual double model(const GModels &models, const GEvent &event, GVector *gradients=NULL) const
Return model value and (optionally) gradients.
const int & rows(void) const
Return number of matrix rows.
const GEnergy & emin(void) const
Return minimum energy of all intervals.
Definition: GEbounds.hpp:187
Abstract data model base class interface definition.
GObservation(void)
Void constructor.
#define G_EVENTS
virtual bool use_event_for_likelihood(const int &index) const
Check whether bin should be used for likelihood analysis.
virtual double likelihood_poisson_unbinned(const GModels &models, GVector *gradients, GMatrixSparse *curvature, double *npred) const
Evaluate log-likelihood function for Poisson statistic and unbinned analysis (version with working ar...
GVector log(const GVector &vector)
Computes natural logarithm of vector elements.
Definition: GVector.cpp:1274
const std::string & name(void) const
Return observation name.
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.
virtual double likelihood_poisson_binned(const GModels &models, GVector *gradients, GMatrixSparse *curvature, double *npred) const
Evaluate log-likelihood function for Poisson statistic and binned analysis (version with working arra...
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.
Abstract response base class definition.
const double & factor_min(void) const
Return parameter minimum factor boundary.
const double & secs(void) const
Return time in seconds in native reference (TT)
Definition: GTime.hpp:156
#define G_LIKELIHOOD
Sky model class interface definition.
Abstract event container class.
Definition: GEvents.hpp:66
virtual GVector column(const int &column) const
Extract column as vector from matrix.
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.
void stack_init(const int &size=0, const int &entries=0)
Initialises matrix filling stack.
Exception handler interface definition.
int npars(void) const
Return number of model parameters in container.
Definition: GModels.cpp:853
const GModel * m_model
Model.
void fixed_iter(const int &iter)
Set fixed number of iterations.
Definition: GIntegral.hpp:183
std::string toupper(const std::string &s)
Convert string to upper case.
Definition: GTools.cpp:941
std::vector< std::string > m_pars_with_gradients
GVector exp(const GVector &vector)
Computes exponential of vector elements.
Definition: GVector.cpp:1232
const int & size(void) const
Return size of vector.
Definition: GVector.hpp:178
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:199
void copy_members(const GObservation &obs)
Copy class members.
Abstract instrument response base class.
Definition: GResponse.hpp:77
virtual int number(void) const =0
bool has_gradient(const GModel &model, const GModelPar &par) const
Check whether a model parameter has an analytical gradient.
Integration class interface definition.
virtual const GResponse * response(void) const =0
Abstract event bin container class.
Definition: GEventCube.hpp:46
const int & columns(void) const
Return number of matrix columns.
int size(void) const
Return number of models in container.
Definition: GModels.hpp:259
const GObservation * m_parent
Observation.
const double & factor_value(void) const
Return parameter factor value.
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:489