Merge branch '1863_CovMat' into integration

# Conflicts: # ChangeLog # NEWS

Merge branch '1863_CovMat' into integration
# Conflicts: # ChangeLog # NEWS
Jürgen Knödlseder
2 parents ad9fa30c 78192ad8
Showing 4 changed files with 82 additions and 2 deletions
ChangeLog
NEWS
include/GObservations.hpp
src/obs/GObservations.cpp
-2016-10-16
+2016-10-17
         * Version 1.2.0 released
         =========================
+        Add GObservations::save_covmat() method (#1863)
         Add GFft and GFftWavetable classes (#1731)
         Add GNdarray class (#1768)
         Increase precision of gammalib::str() floating point functions
 New Features and Important Changes in GammaLib 1.2.0
-16 October 2016
+17 October 2016
 1. Introduction
@@ -54,6 +54,7 @@ The following methods have been added:
 - GTime::now()
 - GFits::publish()
 - gammalib::rstrip_chars()
+- GObservations::save_covmat()
 The following methods have been removed:
 - <all classes>::eval_gradients()
@@ -179,6 +180,9 @@ times that depend on the geographic longitude. Finally, a now() method has
 been added to the GTime class that sets the time to the current time (see
 #1856).
+The GObservations::save_covmat() method computes and saves the covariance
+matrix into a FITS file (#1863).
+
 11. Optimization module
 -----------------------
@@ -129,6 +129,7 @@ public:
     void                optimize(GOptimizer& opt);
     void                errors(GOptimizer& opt);
     void                errors_hessian(void);
+    void                save_covmat(const GFilename& filename) const;
     void                eval(void);
     double              logL(void) const;
     int                 nobserved(void) const;
@@ -35,6 +35,9 @@
 #include "GObservationRegistry.hpp"
 #include "GMatrixSparse.hpp"
 #include "GOptimizer.hpp"
+#include "GFitsBinTable.hpp"
+#include "GFitsTableStringCol.hpp"
+#include "GFitsTableDoubleCol.hpp"
 /* __ Method name definitions ____________________________________________ */
 #define G_AT                                        "GObservations::at(int&)"
@@ -837,6 +840,77 @@ void GObservations::errors_hessian(void)
 /***********************************************************************//**
+ * @brief Save covariance matrix into FITS file
+ *
+ * @param[in] filename FITS filename.
+ *
+ * Saves the covariance matrix into FITS file. The covariance matrix is
+ * calculated as the inverse of the curvature matrix.
+ ***************************************************************************/
+void GObservations::save_covmat(const GFilename& filename) const
+{
+    // Calculate covariance matrix (circumvent const correctness)
+    GMatrixSparse covmat =
+              const_cast<GObservations*>(this)->m_fct.curvature()->invert();
+
+    // Create covariance matrix entry names
+    std::vector<std::string> covmat_entries;
+    for (int i = 0 ; i < m_models.size(); ++i) {
+	    for (int j = 0; j < m_models[i]->size(); ++j) {
+            covmat_entries.push_back(m_models[i]->at(j).name() + "(" +
+                                     m_models[i]->name() + ")");
+        }
+    }
+
+    // Create binary table and columns
+    int size = covmat_entries.size();
+    GFitsBinTable       covmat_table;
+    GFitsTableStringCol par("Parameters", 1, 50, size);
+    GFitsTableDoubleCol cov("Covariance", 1, size*size);
+
+    // Fill tables. We need to take the scale factors into account since the
+    // covariance matrix is for the unscaled parameters
+    int counter = 0;
+    for (int i = 0; i < size; ++i) {
+        par(0, i) = covmat_entries[i];
+        for (int j = 0; j < size; ++j) {
+            cov(0, counter) = covmat(i,j) * m_models.pars()[i]->scale()
+		                                  * m_models.pars()[j]->scale();
+            ++counter;
+        }
+    }
+
+    // Set dimension for covariance matrix column
+    std::vector<int> dim;
+    dim.push_back(size);
+    dim.push_back(size);
+    cov.dim(dim);
+
+    // Append columns to table
+    covmat_table.append(par);
+    covmat_table.append(cov);
+
+    // Set extension name
+    covmat_table.extname("Covariance Matrix");
+
+    // Allocate FITS object
+    GFits fits;
+
+    // Append covariance matrix table to FITS object
+    fits.append(covmat_table);
+
+    // Save FITS file to disk
+    fits.saveto(filename, true);
+
+    // Close FITS object
+    fits.close();
+
+    // Return
+    return;
+}
+
+
+/***********************************************************************//**
  * @brief Evaluate function
  *
  * Evaluates the likelihood funtion at the actual set of parameters.