From d8375b6b897d01ee44ea4ed91d012eb586a85b6c Mon Sep 17 00:00:00 2001
From: Alexander Ziegler <alexander.ziegler@fau.de>
Date: Thu, 6 Oct 2016 17:35:03 +0200
Subject: [PATCH] Add function to write out covariance matrix in FITS file

---
 include/GObservations.hpp |  1 +
 src/obs/GObservations.cpp | 66 ++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
 2 files changed, 67 insertions(+), 0 deletions(-)

diff --git a/include/GObservations.hpp b/include/GObservations.hpp
index dc88904..4578e69 100644
--- a/include/GObservations.hpp
+++ b/include/GObservations.hpp
@@ -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;
diff --git a/src/obs/GObservations.cpp b/src/obs/GObservations.cpp
index 3310f2b..135cebd 100644
--- a/src/obs/GObservations.cpp
+++ b/src/obs/GObservations.cpp
@@ -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
--
libgit2 0.21.4