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