Integrate code

Jürgen Knödlseder
1 parent da257d99
Showing 7 changed files with 61 additions and 39 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-09-26
+2016-10-17
         * Version 1.2.0 released
         ========================
+        Add outcovmat parameter to ctlike (#1863)
         Fix ctbin segfault bug for varying energy thresholds (#1850)
         Make csresmap enumbins parameter hidden (#1844)
         Rename ctskymap::map() to ctskymap::skymap() (#1828)
 New Features and Important Changes in ctools 1.2.0
-26 September 2016
+17 October 2016
 Introduction
@@ -93,6 +93,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
 ------------------------------------------
@@ -49,7 +49,7 @@ 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.
 In addition, the entire output of the covariance matrix to a separate FITS file
-can be queried (parameter ``covmat``). The parameter order given in the 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
@@ -88,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).
@@ -97,11 +100,7 @@ General parameters
 ``(fix_spat_for_ts = no) [boolean]``
     Fix spatial parameters for TS computation.
-``(covmat = NONE) [string]``
-    Output FITS file to store covariance matrix.
-
- 
 Standard parameters
 -------------------
@@ -228,19 +228,13 @@ void ctlike::run(void)
     // Optimize model parameters using LM optimizer
     optimize_lm();
-    // Save covariance matrix if valid filename is passed
-    if (is_valid_filename(m_covmat)) {
-
-        // Write header
-        log_header1(TERSE, "Saving covariance matrix");
-
-        // Log filename
-        log_value(NORMAL, "Covariance matrix file", m_covmat.url());
-
-        // Save covarianve matrix
-        m_obs.save_covmat(m_covmat);
-    }
+    // 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();
@@ -318,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;
 }
@@ -326,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)) {
@@ -344,6 +343,7 @@ void ctlike::save(void)
         // Write results out as XML model
         m_obs.models().save(m_outmodel);
+
     }
     // ... otherwise signal that file was not saved
@@ -351,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;
 }
@@ -405,12 +421,12 @@ void ctlike::get_parameters(void)
     m_apply_edisp     = (*this)["edisp"].boolean();
     m_fix_spat_for_ts = (*this)["fix_spat_for_ts"].boolean();
     m_chatter         = static_cast<GChatter>((*this)["chatter"].integer());
-    m_covmat          = (*this)["covmat"].filename();
     // 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
@@ -543,7 +559,7 @@ void ctlike::init_members(void)
 {
     // Initialise members
     m_outmodel.clear();
-    m_covmat.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
@@ -574,7 +590,7 @@ void ctlike::copy_members(const ctlike&amp; app)
     // Copy attributes
     m_refit           = app.m_refit;
     m_outmodel        = app.m_outmodel;
-    m_covmat          = app.m_covmat;
+    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,10 +70,10 @@ 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
-    GFilename     m_covmat;          //!< Covariance matrix output FITS file name
     // Members
     int           m_max_iter;   //!< Maximum number of iterations
@@ -5,7 +5,7 @@
 # ------------------------------------------------------------------------- #
 #                                                                           #
 #  File:      ctlike.par                                                    #
-#  Version:   1.1.0                                                         #
+#  Version:   1.2.0                                                         #
 #                                                                           #
 #  Author:    Jurgen Knodlseder                                             #
 #             Jurgen.Knodlseder@irap.omp.eu                                 #
@@ -17,17 +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"
-covmat,          f, h, "NONE",,, "Output covariance matrix FITS 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