Rolf Buehler / gammalib

Commit 3e1ac4f11250c90a763d9d2ba865015d46625d3f

Authored by Jürgen Knödlseder
1 parent 59758162

Add missing GObservations::likelihood::covariance() and GObservations::likelihoo… 

…d::save() methods to Python interface and avoid segfault in GObservations::likelihood::covariance() method for non-allocated curvature matrix
Inline Side-by-side
Showing 2 changed files with 23 additions and 11 deletions
pyext/GObservations.i
  View file @3e1ac4f
... ... @@ -94,6 +94,8 @@ public:
94 94 void set(GObservations* obs);
95 95 double npred(void) const;
96 96 GMatrixSparse hessian(const GOptimizerPars& pars);
  97 + GMatrixSparse covariance(void) const;
  98 + void save(const GFilename& filename) const;
97 99 };
98 100 %nestedworkaround GObservations::likelihood;
99 101 %{
... ...
src/obs/GObservations_likelihood.cpp
  View file @3e1ac4f
... ... @@ -388,7 +388,6 @@ GMatrixSparse GObservations::likelihood::hessian(const GOptimizerPars& pars)
388 388 //const double step_tolerance = 0.1;
389 389 //const double gradient_tolerance = 0.02;
390 390  
391   -
392 391 // Create working copy of parameters
393 392 GOptimizerPars wrk_pars = pars;
394 393  
... ... @@ -598,22 +597,33 @@ GMatrixSparse GObservations::likelihood::hessian(const GOptimizerPars& pars)
598 597 * that are stored with the observations and multiplies the covariance
599 598 * matrix with the scaling factors of the model parameters. Hence the
600 599 * covariance matrix is covariance with respect to the true model values.
  600 + *
  601 + * If the curvature matrix does not exist, an empty covariance matrix is
  602 + * returned.
601 603 ***************************************************************************/
602 604 GMatrixSparse GObservations::likelihood::covariance(void) const
603 605 {
604   - // Compute covariance matrix (circumvent const correctness)
605   - GMatrixSparse covmat =
606   - const_cast<GObservations::likelihood*>(this)->curvature()->invert();
  606 + // Initialise covariance matrix
  607 + GMatrixSparse covmat;
  608 +
  609 + // Get pointer on curvature matrix. Only continue if pointer is valid
  610 + GMatrixSparse* curvature = const_cast<GObservations::likelihood*>(this)->curvature();
  611 + if (curvature != NULL) {
  612 +
  613 + // Compute covariance matrix
  614 + covmat = curvature->invert();
607 615  
608   - // Get model parameters
609   - GOptimizerPars pars = m_this->m_models.pars();
  616 + // Get model parameters
  617 + GOptimizerPars pars = m_this->m_models.pars();
610 618  
611   - // Multiply covariance matrix elements with scale factors
612   - for (int row = 0; row < covmat.rows(); ++row) {
613   - for (int col = 0; col < covmat.columns(); ++col) {
614   - covmat(row,col) *= pars[row]->scale() * pars[col]->scale();
  619 + // Multiply covariance matrix elements with scale factors
  620 + for (int row = 0; row < covmat.rows(); ++row) {
  621 + for (int col = 0; col < covmat.columns(); ++col) {
  622 + covmat(row,col) *= pars[row]->scale() * pars[col]->scale();
  623 + }
615 624 }
616   - }
  625 +
  626 + } // endif: curvature matrix was valid
617 627  
618 628 // Return covariance matrix
619 629 return covmat;
... ...