GammaLib  2.0.0
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Macros Pages
GObservations_likelihood.cpp
Go to the documentation of this file.
1 /***************************************************************************
2  * GObservations_likelihood.cpp - Likelihood function class *
3  * ----------------------------------------------------------------------- *
4  * copyright (C) 2009-2021 by Juergen Knoedlseder *
5  * ----------------------------------------------------------------------- *
6  * *
7  * This program is free software: you can redistribute it and/or modify *
8  * it under the terms of the GNU General Public License as published by *
9  * the Free Software Foundation, either version 3 of the License, or *
10  * (at your option) any later version. *
11  * *
12  * This program is distributed in the hope that it will be useful, *
13  * but WITHOUT ANY WARRANTY; without even the implied warranty of *
14  * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the *
15  * GNU General Public License for more details. *
16  * *
17  * You should have received a copy of the GNU General Public License *
18  * along with this program. If not, see <http://www.gnu.org/licenses/>. *
19  * *
20  ***************************************************************************/
21 /**
22  * @file GObservations_likelihood.cpp
23  * @brief Likelihood function class implementation
24  * @author Juergen Knoedlseder
25  */
26 
27 /* __ Includes ___________________________________________________________ */
28 #ifdef HAVE_CONFIG_H
29 #include <config.h>
30 #endif
31 #include "GObservations.hpp"
32 #include "GTools.hpp"
33 #include "GFilename.hpp"
34 #include "GModel.hpp"
35 #include "GEvent.hpp"
36 #include "GEventList.hpp"
37 #include "GEventCube.hpp"
38 #include "GEventBin.hpp"
39 #include "GFitsBinTable.hpp"
40 #include "GFitsTableStringCol.hpp"
41 #include "GFitsTableDoubleCol.hpp"
42 #include "GCsv.hpp"
43 
44 /* __ OpenMP section _____________________________________________________ */
45 #ifdef _OPENMP
46 #include <omp.h>
47 #endif
48 
49 /* __ Method name definitions ____________________________________________ */
50 #define G_EVAL "GObservations::likelihood::eval(GOptimizerPars&)"
51 
52 /* __ Macros _____________________________________________________________ */
53 
54 /* __ Coding definitions _________________________________________________ */
55 //#define G_USE_HESSIAN
56 
57 /* __ Debug definitions __________________________________________________ */
58 //#define G_EVAL_TIMING //!< Perform optimizer timing
59 //#define G_EVAL_DEBUG //!< Perform optimizer debugging
60 //#define G_HESSIAN //!< Debug Hessian computation
61 
62 /* __ Prototypes _________________________________________________________ */
63 
64 
65 /*==========================================================================
66  = =
67  = Constructors/destructors =
68  = =
69  ==========================================================================*/
70 
71 /***********************************************************************//**
72  * @brief Void constructor
73  ***************************************************************************/
75 {
76  // Initialise members
77  init_members();
78 
79  // Return
80  return;
81 }
82 
83 
84 /***********************************************************************//**
85  * @brief Observations constructor
86  *
87  * @param[in] obs Observations container pointer.
88  *
89  * Constructs optimizer from GObservations container. The method copies the
90  * pointer to the observation container in the m_this member, making the
91  * observation container accessible to the optimizer class.
92  ***************************************************************************/
94 {
95  // Initialise members
96  init_members();
97 
98  // Set object
99  m_this = obs;
100 
101  // Return
102  return;
103 }
104 
105 
106 
107 /***********************************************************************//**
108  * @brief Copy constructor
109  *
110  * @param[in] fct Optimizer function.
111  ***************************************************************************/
113  GOptimizerFunction(fct)
114 {
115  // Initialise members
116  init_members();
117 
118  // Copy members
119  copy_members(fct);
120 
121  // Return
122  return;
123 }
124 
125 
126 /***********************************************************************//**
127  * @brief Destructor
128  ***************************************************************************/
130 {
131  // Free members
132  free_members();
133 
134  // Return
135  return;
136 }
137 
138 
139 /*==========================================================================
140  = =
141  = Operators =
142  = =
143  ==========================================================================*/
144 
145 /***********************************************************************//**
146  * @brief Assignment operator
147  *
148  * @param[in] fct Likelihood function.
149  * @return Likelihood function.
150  ***************************************************************************/
152 {
153  // Execute only if object is not identical
154  if (this != &fct) {
155 
156  // Copy base class members
158 
159  // Free members
160  free_members();
161 
162  // Initialise private members
163  init_members();
164 
165  // Copy members
166  copy_members(fct);
167 
168  } // endif: object was not identical
169 
170  // Return this object
171  return *this;
172 }
173 
174 
175 /*==========================================================================
176  = =
177  = Public methods =
178  = =
179  ==========================================================================*/
180 
181 /***********************************************************************//**
182  * @brief Evaluate log-likelihood function
183  *
184  * @param[in] pars Optimizer parameters.
185  *
186  * @exception GException::invalid_statistics
187  * Invalid optimization statistics encountered.
188  *
189  * This method evaluates the -(log-likelihood) function for parameter
190  * optimization. It handles both binned and unbinned data and supportes
191  * Poisson and Gaussian statistics.
192  * Note that different statistics and different analysis methods
193  * (binned/unbinned) may be combined.
194  ***************************************************************************/
196 {
197  // Timing measurement
198  #if defined(G_EVAL_TIMING)
199  #ifdef _OPENMP
200  double t_start = omp_get_wtime();
201  #else
202  clock_t t_start = clock();
203  #endif
204  #endif
205 
206  // Single loop for common exit point
207  do {
208 
209  // Get number of parameters for allocation of vectors and matrices
210  int npars = pars.size();
211 
212  // Free old memory
213  if (m_gradient != NULL) delete m_gradient;
214  if (m_curvature != NULL) delete m_curvature;
215 
216  // Initialise value, gradient vector and curvature matrix
217  m_value = 0.0;
218  m_npred = 0.0;
219  m_gradient = new GVector(npars);
220  m_curvature = new GMatrixSparse(npars,npars);
221 
222  // Set stack size and number of entries
223  int stack_size = (2*npars > 100000) ? 2*npars : 100000;
224  int max_entries = 2*npars;
225  m_curvature->stack_init(stack_size, max_entries);
226 
227  // Allocate vectors to save working variables of each thread
228  std::vector<GVector*> vect_cpy_grad;
229  std::vector<GMatrixSparse*> vect_cpy_curvature;
230  std::vector<double*> vect_cpy_value;
231  std::vector<double*> vect_cpy_npred;
232 
233  // Here OpenMP will paralellize the execution. The following code will
234  // be executed by the differents threads. In order to avoid protecting
235  // attributes ( m_value,m_npred, m_gradient and m_curvature), each thread
236  // works with its own working variables (cpy_*). When a thread starts,
237  // we add working variables in a vector (vect_cpy_*). When computation
238  // is finished we just add all elements contain in the vector to the
239  // attributes value.
240  #pragma omp parallel
241  {
242  // Allocate and initialize variable copies for multi-threading
243  GModels cpy_model(m_this->models());
244  GVector* cpy_gradient = new GVector(npars);
245  GMatrixSparse* cpy_curvature = new GMatrixSparse(npars,npars);
246  double* cpy_npred = new double(0.0);
247  double* cpy_value = new double(0.0);
248 
249  // Set stack size and number of entries
250  cpy_curvature->stack_init(stack_size, max_entries);
251 
252  // Push variable copies into vector. This is a critical zone to
253  // avoid multiple thread pushing simultaneously.
254  #pragma omp critical(GObservations_likelihood_eval)
255  {
256  vect_cpy_grad.push_back(cpy_gradient);
257  vect_cpy_curvature.push_back(cpy_curvature);
258  vect_cpy_value.push_back(cpy_value);
259  vect_cpy_npred.push_back(cpy_npred);
260  }
261 
262  // Loop over all observations. The omp for directive will deal
263  // with the iterations on the differents threads.
264  #pragma omp for
265  for (int i = 0; i < m_this->size(); ++i) {
266 
267  // Compute likelihood
268  *cpy_value += m_this->m_obs[i]->likelihood(cpy_model,
269  cpy_gradient,
270  cpy_curvature,
271  cpy_npred);
272 
273  } // endfor: looped over observations
274 
275  // Release stack
276  cpy_curvature->stack_destroy();
277 
278  } // end pragma omp parallel
279 
280  // Now the computation is finished, update attributes.
281  // For each omp section, a thread will be created.
282  #pragma omp parallel sections
283  {
284  #pragma omp section
285  {
286  for (int i = 0; i < vect_cpy_curvature.size() ; ++i) {
287  *m_curvature += *(vect_cpy_curvature.at(i));
288  delete vect_cpy_curvature.at(i);
289  }
290  }
291 
292  #pragma omp section
293  {
294  for (int i = 0; i < vect_cpy_grad.size(); ++i){
295  *m_gradient += *(vect_cpy_grad.at(i));
296  delete vect_cpy_grad.at(i);
297  }
298  }
299 
300  #pragma omp section
301  {
302  for(int i = 0; i < vect_cpy_npred.size(); ++i){
303  m_npred += *(vect_cpy_npred.at(i));
304  delete vect_cpy_npred.at(i);
305  }
306  }
307 
308  #pragma omp section
309  {
310  for (int i = 0; i < vect_cpy_value.size(); ++i){
311  m_value += *(vect_cpy_value.at(i));
312  delete vect_cpy_value.at(i);
313  }
314  }
315  } // end of pragma omp sections
316 
317  // Release stack
318  m_curvature->stack_destroy();
319 
320  } while(0); // endwhile: main loop
321 
322  // Copy over the parameter gradients for all parameters that are
323  // free (so that we can access the gradients from outside)
324  for (int ipar = 0; ipar < pars.size(); ++ipar) {
325  if (pars[ipar]->is_free()) {
326  GOptimizerPar* par = const_cast<GOptimizerPar*>(pars[ipar]);
327  par->factor_gradient((*m_gradient)[ipar]);
328  }
329  }
330 
331  // Optionally use Hessian instead of curvature matrix
332  #if defined(G_USE_HESSIAN)
333  *m_curvature = hessian(pars);
334  #endif
335 
336  // Optionally dump gradient and curvature matrix
337  #if defined(G_EVAL_DEBUG)
338  std::cout << *m_gradient << std::endl;
339  for (int i = 0; i < pars.size(); ++i) {
340  for (int j = 0; j < pars.size(); ++j) {
341  std::cout << (*m_curvature)(i,j) << " ";
342  }
343  std::cout << std::endl;
344  }
345  #endif
346 
347  // Timing measurement
348  #if defined(G_EVAL_TIMING)
349  #ifdef _OPENMP
350  double t_elapse = omp_get_wtime()-t_start;
351  #else
352  double t_elapse = (double)(clock() - t_start) / (double)CLOCKS_PER_SEC;
353  #endif
354  std::cout << "GObservations::optimizer::eval: CPU usage = "
355  << t_elapse << " sec" << std::endl;
356  #endif
357 
358  // Return
359  return;
360 }
361 
362 
363 /***********************************************************************//**
364  * @brief Compute Hessian matrix
365  *
366  * @param[in] pars Optimizer parameters.
367  *
368  * @return Hessian matrix.
369  ***************************************************************************/
371 {
372  // Set strategy constants (low)
373  //const int ncyles = 3;
374  //const double step_tolerance = 0.5;
375  //const double gradient_tolerance = 0.1;
376 
377  // Set strategy constants (medium)
378  const int ncyles = 5;
379  const double step_tolerance = 0.3;
380  const double gradient_tolerance = 0.05;
381 
382  // Set strategy constants (high)
383  //const int ncyles = 7;
384  //const double step_tolerance = 0.1;
385  //const double gradient_tolerance = 0.02;
386 
387  // Create working copy of parameters
388  GOptimizerPars wrk_pars = pars;
389 
390  // Get number of parameters
391  int npars = wrk_pars.size();
392 
393  // Allocate Hessian matrix
394  GMatrixSparse hessian(npars, npars);
395 
396  // Find out machine precision
397  double eps = 0.1;
398  while (1.0+eps != 1.0) {
399  eps *= 0.5;
400  }
401  double eps2 = 2.0 * std::sqrt(eps);
402 
403  // Function value
404  eval(wrk_pars);
405  double f = value();
406 
407  // Compute aimsag
408  double aimsag = std::sqrt(eps2) * std::abs(f);
409 
410  // Diagonal elements
411  std::vector<double> g2(npars, 0.0);
412  std::vector<double> grd(npars, 0.0);
413  std::vector<double> dir(npars, 0.0);
414  std::vector<double> yy(npars, 0.0);
415 
416  // Loop over parameters
417  for (int i = 0; i < npars; ++i) {
418 
419  // Get parameter
420  GOptimizerPar* par = wrk_pars[i];
421 
422  // Interrupt if parameter is fixed
423  if (par->is_fixed()) {
424  hessian(i,i) = 0.0;
425  continue;
426  }
427 
428  // Setup step size
429  double xtf = par->factor_value();
430  double dmin = 8.0 * eps2 * std::abs(xtf);
431  double d = 0.000001;
432  if (d < dmin) {
433  d = dmin;
434  }
435 
436  // Loop over cycles
437  for (int icyc = 0; icyc < ncyles; ++icyc) {
438  //for (int icyc = 0; icyc < 1; ++icyc) {
439 
440  // Initialise
441  double sag = 0.0;
442  double fs1 = 0.0; // right-hand side
443  double fs2 = 0.0; // left-hand side
444 
445  // Compute gradient
446  for (int multpy = 0; multpy < 5; ++multpy) {
447  //for (int multpy = 0; multpy < 1; ++multpy) {
448 
449  // Compute right-hand side
450  par->factor_value(xtf + d);
451  eval(wrk_pars);
452  fs1 = value();
453 
454  // Compute left-hand side
455  par->factor_value(xtf - d);
456  eval(wrk_pars);
457  fs2 = value();
458 
459  // Recover current value
460  par->factor_value(xtf);
461 
462  // Compute sag
463  sag = 0.5 * (fs1 + fs2 - 2.0*f);
464 
465  // Break if sag is okay
466  if (std::abs(sag) > eps2 || sag == 0.0) {
467  break;
468  }
469 
470  // ... otherwise increase step size
471  d *= 10.0;
472 
473  } // endfor
474 
475  // Save old step size and second derivative
476  double dlast = d;
477  double g2bfor = g2[i];
478 
479  // Compute parameter derivatives and store step size and
480  // function value
481  g2[i] = 2.0 * sag/(d*d);
482  grd[i] = (fs1-fs2)/(2.*d);
483  dir[i] = d;
484  yy[i] = fs1;
485 
486  // Compute a new step size based on the aimed sag
487  if (sag != 0.0) {
488  d = std::sqrt(2.0*aimsag/std::abs(g2[i]));
489  }
490  if (d < dmin) {
491  d = dmin;
492  }
493  /*
494  else if (par->factor_value()+d > par->factor_max()) {
495  d = dmin;
496  }
497  else if (par->factor_value()-d > par->factor_min()) {
498  d = dmin;
499  }
500  */
501 
502  // Check if converged
503  if (std::abs((d-dlast)/d) < step_tolerance) {
504  break;
505  }
506  if (std::abs((g2[i]-g2bfor)/g2[i]) < gradient_tolerance) {
507  break;
508  }
509  d = std::min(d, 10.*dlast);
510  d = std::max(d, 0.1*dlast);
511 
512  } // endfor: cycles
513 
514  // Set diagonal element
515  hessian(i,i) = g2[i];
516 
517  } // endfor: looped over all parameters
518 
519  // Debug dump
520  #if defined(G_HESSIAN)
521  std::cout << "GObservations::likelihood::hessian: ";
522  std::cout << "deltas and gradients:" << std::endl;
523  for (int i = 0; i < npars; ++i) {
524  std::cout << dir[i] << " ";
525  std::cout << grd[i] << std::endl;
526  }
527  #endif
528 
529  // Compute off-diagonal elements
530  for (int i = 0; i < npars; ++i) {
531 
532  // Get parameter 1
533  GOptimizerPar* par1 = wrk_pars[i];
534  double x1 = par1->factor_value();
535 
536  // Increment parameter 1
537  par1->factor_value(x1 + dir[i]);
538 
539  // Loop over columns
540  for (int j = i+1; j < npars; ++j) {
541 
542  // Get parameter 2
543  GOptimizerPar* par2 = wrk_pars[j];
544  double x2 = par2->factor_value();
545 
546  // Interrupt if parameter is fixed
547  if (par1->is_fixed() || par2->is_fixed()) {
548  hessian(i,j) = 0.0;
549  hessian(j,i) = 0.0;
550  continue;
551  }
552 
553  // Increment parameter 2
554  par2->factor_value(x2 + dir[j]);
555 
556  // Evaluate Hessian element
557  eval(wrk_pars);
558  double fs1 = value();
559  double element = (fs1 + f - yy[i] - yy[j])/(dir[i]*dir[j]);
560 
561  // Store Hessian element
562  hessian(i,j) = element;
563  hessian(j,i) = element;
564 
565  // Restore parameter 2
566  par2->factor_value(x2);
567 
568  } // endfor: looped over columns
569 
570  // Restore parameter 1
571  par1->factor_value(x1);
572 
573  } // endfor: looped over parameters
574 
575  // Debug dump
576  #if defined(G_HESSIAN)
577  std::cout << "GObservations::likelihood::hessian: " << std::endl;
578  std::cout << hessian << std::endl;
579  #endif
580 
581  // Return Hessian
582  return hessian;
583 }
584 
585 
586 /***********************************************************************//**
587  * @brief Compute covariance matrix
588  *
589  * @return Covariance matrix.
590  *
591  * The covariance matrix is calculated as the inverse of the curvature
592  * matrix. The method retrieves the model scaling factors from the models
593  * that are stored with the observations and multiplies the covariance
594  * matrix with the scaling factors of the model parameters. Hence the
595  * covariance matrix is covariance with respect to the true model values.
596  *
597  * If the curvature matrix does not exist, an empty covariance matrix is
598  * returned.
599  ***************************************************************************/
601 {
602  // Initialise covariance matrix
603  GMatrixSparse covmat;
604 
605  // Get pointer on curvature matrix. Only continue if pointer is valid
606  GMatrixSparse* curvature = const_cast<GObservations::likelihood*>(this)->curvature();
607  if (curvature != NULL) {
608 
609  // Compute covariance matrix
610  covmat = curvature->invert();
611 
612  // Get model parameters
613  GOptimizerPars pars = m_this->m_models.pars();
614 
615  // Multiply covariance matrix elements with scale factors
616  for (int row = 0; row < covmat.rows(); ++row) {
617  for (int col = 0; col < covmat.columns(); ++col) {
618  covmat(row,col) *= pars[row]->scale() * pars[col]->scale();
619  }
620  }
621 
622  } // endif: curvature matrix was valid
623 
624  // Return covariance matrix
625  return covmat;
626 }
627 
628 
629 /***********************************************************************//**
630  * @brief Save likelihood fit results into a CSV or FITS file.
631  *
632  * @param[in] filename CSV or FITS filename.
633  *
634  * Saves the likelihood fit results into a CSV or FITS file. The result
635  * format depends on the filename extension. If the extension is `.fits` or
636  * `.fit` the file is written into a FITS file, otherwise it is written into
637  * a CSV file.
638  ***************************************************************************/
639 void GObservations::likelihood::save(const GFilename& filename) const
640 {
641  // Get file type
642  std::string filetype = filename.type();
643 
644  // If file type is a FITS file then write save covariance matrix into a
645  // FITS file, otherwise save the covariance matrix into a CSV file
646  if (filetype == "fits") {
647  save_fits(filename);
648  }
649  else {
650  save_csv(filename);
651  }
652 
653  // Return
654  return;
655 }
656 
657 
658 /*==========================================================================
659  = =
660  = Private methods =
661  = =
662  ==========================================================================*/
663 
664 /***********************************************************************//**
665  * @brief Initialise class members
666  ***************************************************************************/
668 {
669  // Initialise members
670  m_value = 0.0;
671  m_npred = 0.0;
672  m_this = NULL;
673  m_gradient = NULL;
674  m_curvature = NULL;
675 
676  // Return
677  return;
678 }
679 
680 
681 /***********************************************************************//**
682  * @brief Copy class members
683  *
684  * @param[in] fct Optimizer.
685  ***************************************************************************/
687 {
688  // Copy attributes
689  m_value = fct.m_value;
690  m_npred = fct.m_npred;
691  m_this = fct.m_this;
692 
693  // Clone gradient if it exists
694  if (fct.m_gradient != NULL) {
695  m_gradient = new GVector(*fct.m_gradient);
696  }
697 
698  // Clone curvature matrix if it exists
699  if (fct.m_curvature != NULL) {
700  m_curvature = new GMatrixSparse(*fct.m_curvature);
701  }
702 
703  // Return
704  return;
705 }
706 
707 
708 /***********************************************************************//**
709  * @brief Delete class members
710  ***************************************************************************/
712 {
713  // Free members
714  if (m_gradient != NULL) delete m_gradient;
715  if (m_curvature != NULL) delete m_curvature;
716 
717  // Signal free pointers
718  m_gradient = NULL;
719  m_curvature = NULL;
720 
721  // Return
722  return;
723 }
724 
725 
726 /***********************************************************************//**
727  * @brief Save likelihood fit results into CSV file.
728  *
729  * @param[in] filename CSV filename.
730  *
731  * Saves the likelihood fit results into a CSV file. For the moment the
732  * method only writes the covariance matrix.
733  ***************************************************************************/
735 {
736  // Get covariance matrix
737  GMatrixSparse covmat = covariance();
738 
739  // Get covariance matrix row and column names
740  std::vector<std::string> names = covariance_names();
741 
742  // Create binary table and columns
743  int size = names.size();
744 
745  // Create binary table and columns
746  GCsv table(size + 1, size);
747 
748  // Fill CSV header
749  for (int i = 0; i < size; ++i) {
750  table.string(0, i, names[i]);
751  }
752 
753  // Fill CSV table
754  for (int i = 0; i < size; ++i) {
755  for (int j = 0; j < size; j ++) {
756  table.real(i+1, j, covmat(i,j));
757  }
758  }
759 
760  // Save CSV table
761  table.save(filename, " ", true);
762 
763  // Return
764  return;
765 }
766 
767 
768 /***********************************************************************//**
769  * @brief Save likelihood fit results into FITS file.
770  *
771  * @param[in] filename FITS filename.
772  *
773  * Saves the likelihood fit results into a FITS file. For the moment the
774  * method only writes the covariance matrix.
775  ***************************************************************************/
777 {
778  // Get covariance matrix
779  GMatrixSparse covmat = covariance();
780 
781  // Get covariance matrix row and column names
782  std::vector<std::string> names = covariance_names();
783 
784  // Create binary table and columns
785  int size = names.size();
786 
787  // Create binary table and columns
788  GFitsBinTable covmat_table;
789  GFitsTableStringCol par("Parameters", 1, 50, size);
790  GFitsTableDoubleCol cov("Covariance", 1, size*size);
791 
792  // Fill tables
793  int counter = 0;
794  for (int i = 0; i < size; ++i) {
795  par(0, i) = names[i];
796  for (int j = 0; j < size; ++j) {
797  cov(0, counter) = covmat(i,j);
798  ++counter;
799  }
800  }
801 
802  // Set dimension for covariance matrix column
803  std::vector<int> dim;
804  dim.push_back(size);
805  dim.push_back(size);
806  cov.dim(dim);
807 
808  // Append columns to table
809  covmat_table.append(par);
810  covmat_table.append(cov);
811 
812  // Set extension name
813  covmat_table.extname("Covariance Matrix");
814 
815  // Allocate FITS object
816  GFits fits;
817 
818  // Append covariance matrix table to FITS object
819  fits.append(covmat_table);
820 
821  // Save FITS file to disk
822  fits.saveto(filename, true);
823 
824  // Close FITS object
825  fits.close();
826 
827  // Return
828  return;
829 }
830 
831 
832 /***********************************************************************//**
833  * @brief Return covariance matrix row and column names
834  *
835  * @return Covariance matrix row and column names.
836  *
837  * Returns the row and column names of the covariance matrix.
838  ***************************************************************************/
839 std::vector<std::string> GObservations::likelihood::covariance_names(void) const
840 {
841  // Initialise covariance matrix row and column names
842  std::vector<std::string> names;
843 
844  // Create covariance matrix row and column names
845  for (int i = 0 ; i < m_this->m_models.size(); ++i) {
846  for (int j = 0; j < m_this->m_models[i]->size(); ++j) {
847  names.push_back(m_this->m_models[i]->at(j).name() + "(" +
848  m_this->m_models[i]->name() + ")");
849  }
850  }
851 
852  // Return covariance matrix row and column names
853  return names;
854 }
GObservations * m_this
Pointer to GObservations object.
likelihood(void)
Void constructor.
void save(const GFilename &filename, const std::string &sep=" ", const bool &clobber=false) const
Save CSV table.
Definition: GCsv.cpp:519
FITS table double column class interface definition.
double m_value
Function value.
Abstract model base class interface definition.
Optimizer function abstract base class.
const double & factor_gradient(void) const
Return parameter factor gradient.
void stack_destroy(void)
Destroy matrix stack.
Sparse matrix class interface definition.
void save_csv(const GFilename &filename) const
Save likelihood fit results into CSV file.
GVector abs(const GVector &vector)
Computes absolute of vector elements.
Definition: GVector.cpp:1253
Abstract event bin base class definition.
Abstract event bin container class interface definition.
Optimizer parameter container class.
void save(const GFilename &filename) const
Save likelihood fit results into a CSV or FITS file.
GVector * m_gradient
Pointer to gradient vector.
GFitsTableCol * append(const GFitsTableCol &column)
Append column to the table.
Definition: GFitsTable.hpp:147
GMatrixSparse hessian(const GOptimizerPars &pars)
Compute Hessian matrix.
Comma-separated values table class.
Definition: GCsv.hpp:57
GMatrixSparse * m_curvature
Pointer to curvature matrix.
Gammalib tools definition.
GMatrixSparse invert(void) const
Return inverted matrix.
FITS file class.
Definition: GFits.hpp:63
double m_npred
Total number of predicted events.
double min(const GVector &vector)
Computes minimum vector element.
Definition: GVector.cpp:886
void save_fits(const GFilename &filename) const
Save likelihood fit results into FITS file.
Abstract event atom container class interface definition.
void free_members(void)
Delete class members.
FITS table string column.
std::string string(const int &row, const int &col) const
Get string value.
Definition: GCsv.cpp:309
void saveto(const GFilename &filename, const bool &clobber=false)
Saves to specified FITS file.
Definition: GFits.cpp:1293
void copy_members(const GObservations &obs)
Copy class members.
GVector sqrt(const GVector &vector)
Computes square root of vector elements.
Definition: GVector.cpp:1358
Model container class.
Definition: GModels.hpp:152
Filename class.
Definition: GFilename.hpp:62
const int & rows(void) const
Return number of matrix rows.
Abstract event base class definition.
std::vector< std::string > covariance_names(void) const
Return covariance matrix row and column names.
FITS table string column class interface definition.
void init_members(void)
Initialise class members.
likelihood & operator=(const likelihood &fct)
Assignment operator.
bool is_fixed(void) const
Signal if parameter is fixed.
void dim(const std::vector< int > &dim)
Set column dimension.
int size(void) const
Return number of observations in container.
const std::string & extname(void) const
Return extension name.
Definition: GFitsHDU.hpp:162
double max(const GVector &vector)
Computes maximum vector element.
Definition: GVector.cpp:915
void eval(void)
Evaluate function.
Observation container class.
std::string type(void) const
Return file type.
Definition: GFilename.cpp:264
virtual void eval(const GOptimizerPars &pars)
Evaluate log-likelihood function.
double real(const int &row, const int &col) const
Get double precision value.
Definition: GCsv.cpp:324
GMatrixSparse covariance(void) const
Compute covariance matrix.
FITS binary table class.
void stack_init(const int &size=0, const int &entries=0)
Initialises matrix filling stack.
void copy_members(const likelihood &fct)
Copy class members.
GFitsHDU * append(const GFitsHDU &hdu)
Append HDU to FITS file.
Definition: GFits.cpp:678
FITS binary table class definition.
Vector class.
Definition: GVector.hpp:46
int size(void) const
Return number of parameters in container.
void free_members(void)
Delete class members.
void init_members(void)
Initialise class members.
Observations container class interface definition.
const int & columns(void) const
Return number of matrix columns.
void close(void)
Close FITS file.
Definition: GFits.cpp:1342
FITS table double column.
Filename class interface definition.
Comma-separated values table class definition.
const double & factor_value(void) const
Return parameter factor value.
virtual GOptimizerFunction & operator=(const GOptimizerFunction &fct)
Assignment operator.
Optimizer parameter class.