Merge branch '1863_CovMat' into integration

# Conflicts: # ChangeLog # NEWS

Merge branch '1863_CovMat' into integration
# Conflicts: # ChangeLog # NEWS
Jürgen Knödlseder
2 parents 2ccf1353 6cea5977
Showing 7 changed files with 68 additions and 22 deletions
AUTHORS
ChangeLog
NEWS
doc/source/users/reference_manual/ctlike.rst
src/ctlike/ctlike.cpp
src/ctlike/ctlike.hpp
src/ctlike/ctlike.par
@@ -13,3 +13,5 @@ Michael Mayer (michael.mayer@physik.hu-berlin.de)
 Stefan Ohm (stefan.ohm81@gmail.com)
 David Sanchez (david.sanchez@lapp.in2p3.fr)
 Hubert Siejkowski (h.siejkowski@cyfronet.pl)
+Alexander Ziegler (alexander.ziegler@fau.de)
+
-2016-10-07
+2016-10-17
         * Version 1.2.0 released
         ========================
+        Add outcovmat parameter to ctlike (#1863)
         Add publish parameter and methods to ctobssim, ctselect, cslightcrv and csspec (#1871)
         Add csviscube script (#1853)
         Add zenith and azimuth angle information to IRF caldb (#1853)
 New Features and Important Changes in ctools 1.2.0
-7 October 2016
+17 October 2016
 Introduction
@@ -99,6 +99,10 @@ of the tool (#1835).
 The class now derives from ctlikelihood (#1828).
+The hidden parameter "outcovmat" has been added that by default is set to "NONE".
+If instead a FITS filename is given, ctlike will write the coraviance matrix of
+the fit into a FITS file (#1863).
+
 cttsmap - Generation of Test Statistic map
 ------------------------------------------
@@ -48,9 +48,13 @@ the spatial model parameters can be kept fixed during that computation
 ctlike generates an output model XML file that contains the values of the 
 best fitting model parameters. For all free parameters, an ``error`` attribute
 is added that provides the statistical uncertainty in the parameter estimate.
-If computation of the Test Statistics has been requested for a model component,
-a ``ts`` attribute providing the Test Statistics value is added. The output
-model can be used as an input model for other ctools.
+In addition, the entire output of the covariance matrix to a separate FITS file
+can be queried (parameter ``outcovmat``). The parameter order given in the file
+corresponds to the arrangement of the covariance matrix entries. If computation
+of the Test Statistics has been requested for a model component, a ``ts``
+attribute providing the Test Statistics value is added. The output model can be
+used as an input model for other ctools.
+
 General parameters
@@ -84,6 +88,9 @@ General parameters
     Output model XML file with values and uncertainties updated by
     the maximum likelihood fit.
+``(outcovmat = NONE) [string]``
+    Output FITS file to store covariance matrix.
+
 ``(stat = POISSON) [string]``
     Fitting statistics (POISSON or GAUSSIAN; only affects binned analysis).
@@ -92,8 +99,7 @@ General parameters
 ``(fix_spat_for_ts = no) [boolean]``
     Fix spatial parameters for TS computation.
- 	 	 
- 	 	 
+
 Standard parameters
 -------------------
@@ -228,6 +228,13 @@ void ctlike::run(void)
     // Optimize model parameters using LM optimizer
     optimize_lm();
+    // Get pointer to curvature matrix
+    GMatrixSparse* curvature_ptr =
+        const_cast<GObservations::likelihood&>(m_obs.function()).curvature();
+
+    // Store copy of curvature matrix
+    GMatrixSparse curvature = *curvature_ptr;
+    
     // Store Npred
     double npred = m_obs.npred();
@@ -305,6 +312,9 @@ void ctlike::run(void)
     // Restore energy dispersion flags of all CTA observations
     restore_edisp(m_obs, save_edisp);
+    // Restore curvature matrix
+    *curvature_ptr = curvature;
+
     // Return
     return;
 }
@@ -313,15 +323,17 @@ void ctlike::run(void)
 /***********************************************************************//**
  * @brief Save results
  *
- * This method saves the fit results in a XML file.
+ * This method saves the fit results into a XML file and a FITS file. If
+ * the filename parameters are "NONE", the files are not saved.
  ***************************************************************************/
 void ctlike::save(void)
 {
     // Write header
     log_header1(TERSE, "Save results");
-    // Get output filename
-    m_outmodel = (*this)["outmodel"].filename();
+    // Get output filenames
+    m_outmodel  = (*this)["outmodel"].filename();
+    m_outcovmat = (*this)["outcovmat"].filename();
     // Save only if filename is valid
     if (is_valid_filename(m_outmodel)) {
@@ -339,6 +351,22 @@ void ctlike::save(void)
         log_value(NORMAL, "Model definition file", "NONE");
     }
+    // Save covariance matrix if filename is valid
+    if (is_valid_filename(m_outcovmat)) {
+
+        // Log filename
+        log_value(NORMAL, "Covariance matrix file", m_outcovmat.url());
+
+        // Save covarianve matrix
+        m_obs.function().save(m_outcovmat);
+
+    }
+
+    // ... otherwise signal that no covariance matrix was not saved
+    else {
+        log_value(NORMAL, "Covariance matrix file", "NONE");
+    }
+
     // Return
     return;
 }
@@ -397,7 +425,8 @@ void ctlike::get_parameters(void)
     // Optionally read ahead parameters so that they get correctly
     // dumped into the log file
     if (read_ahead()) {
-        m_outmodel = (*this)["outmodel"].filename();
+        m_outmodel  = (*this)["outmodel"].filename();
+        m_outcovmat = (*this)["outcovmat"].filename();
     }
     // Set optimizer logger
@@ -530,6 +559,7 @@ void ctlike::init_members(void)
 {
     // Initialise members
     m_outmodel.clear();
+    m_outcovmat.clear();
     m_refit           = false;
     m_max_iter        = 100;   // Set maximum number of iterations
     m_max_stall       = 10;    // Set maximum number of stalls
@@ -560,6 +590,7 @@ void ctlike::copy_members(const ctlike&amp; app)
     // Copy attributes
     m_refit           = app.m_refit;
     m_outmodel        = app.m_outmodel;
+    m_outcovmat       = app.m_outcovmat;
     m_max_iter        = app.m_max_iter;
     m_max_stall       = app.m_max_stall;
     m_logL            = app.m_logL;
@@ -70,6 +70,7 @@ protected:
     // User parameters
     bool          m_refit;           //!< Refitting
     GFilename     m_outmodel;        //!< Source model output XML file name
+    GFilename     m_outcovmat;       //!< Covariance matrix output FITS file name
     bool          m_apply_edisp;     //!< Apply energy dispersion?
     bool          m_fix_spat_for_ts; //!< Fix spatial parameters for TS computation?
     GChatter      m_chatter;         //!< Chattiness
@@ -5,7 +5,7 @@
 # ------------------------------------------------------------------------- #
 #                                                                           #
 #  File:      ctlike.par                                                    #
-#  Version:   1.1.0                                                         #
+#  Version:   1.2.0                                                         #
 #                                                                           #
 #  Author:    Jurgen Knodlseder                                             #
 #             Jurgen.Knodlseder@irap.omp.eu                                 #
@@ -17,16 +17,17 @@
 #
 # General parameters
 #===================
-inobs,           f, a, "events.fits",,, "Input event list, counts cube or observation definition XML file"
-inmodel,         f, a, "$CTOOLS/share/models/crab.xml",,, "Input model XML file"
-expcube,         f, a, "NONE",,, "Input exposure cube file (only needed for stacked analysis)"
-psfcube,         f, a, "NONE",,, "Input PSF cube file (only needed for stacked analysis)"
-edispcube,       f, a, "NONE",,, "Input energy dispersion cube file (only needed for stacked analysis)"
-bkgcube,         f, a, "NONE",,, "Input background cube file (only needed for stacked analysis)"
-caldb,           s, a, "prod2",,, "Calibration database"
-irf,             s, a, "South_0.5h",,, "Instrument response function"
-edisp,           b, h, no,,, "Apply energy dispersion?"
-outmodel,        f, a, "crab_results.xml",,, "Output model XML file"
+inobs,     f, a, "events.fits",,, "Input event list, counts cube or observation definition XML file"
+inmodel,   f, a, "$CTOOLS/share/models/crab.xml",,, "Input model XML file"
+expcube,   f, a, "NONE",,, "Input exposure cube file (only needed for stacked analysis)"
+psfcube,   f, a, "NONE",,, "Input PSF cube file (only needed for stacked analysis)"
+edispcube, f, a, "NONE",,, "Input energy dispersion cube file (only needed for stacked analysis)"
+bkgcube,   f, a, "NONE",,, "Input background cube file (only needed for stacked analysis)"
+caldb,     s, a, "prod2",,, "Calibration database"
+irf,       s, a, "South_0.5h",,, "Instrument response function"
+edisp,     b, h, no,,, "Apply energy dispersion?"
+outmodel,  f, a, "crab_results.xml",,, "Output model XML file"
+outcovmat, f, h, "NONE",,, "Output covariance matrix FITS file"
 #
 # Maximisation parameters