Add function to write out covariance matrix in FITS file

Alexander Ziegler
1 parent 6f53a5ab
Showing 2 changed files with 67 additions and 0 deletions
include/GObservations.hpp
src/obs/GObservations.cpp
@@ -129,6 +129,7 @@ public:
     void                optimize(GOptimizer& opt);
     void                errors(GOptimizer& opt);
     void                errors_hessian(void);
+    void                save_covmat(const GFilename& filename);
     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&)"
@@ -835,6 +838,69 @@ void GObservations::errors_hessian(void)
     return;
 }
  
+/***********************************************************************//**
+ * @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) 
+{
+    // Calculate covariance matrix
+    GMatrixSparse covmat = 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() + ")" );
+	}
+    }
+
+    // Allocate FITS object
+    GFits fits;
+
+    // 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
+    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);
+            ++counter;
+        }
+    }
+
+    // Set dimension
+    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);
+
+    // Append table to FITS object
+    fits.append(covmat_table);
+
+    // Save FITS file to disk
+    fits.saveto(filename, true);
+
+    // Close FITS object
+    fits.close();
+}
+
  
 /***********************************************************************//**
  * @brief Evaluate function