Improved cstsmapmerge unit testing and updated the test cases, started to adapt …

…cttsmap and cstsmapmerge to the new logging conventions

Improved cstsmapmerge unit testing and updated the test cases, started to adapt …
…cttsmap and cstsmapmerge to the new logging conventions
Michael Mayer
1 parent bc55f554
Showing 7 changed files with 102 additions and 98 deletions
cscripts/cstsmapmerge.py
doc/source/users/reference_manual/cttsmap.rst
src/cttsmap/cttsmap.cpp
src/cttsmap/cttsmap.par
test/data/tsmap_0.fits
test/data/tsmap_1.fits
test/test_cstsmapmerge.py
@@ -46,7 +46,7 @@ class cstsmapmerge(ctools.cscript):
         self._files        = None
         self._in_filename  = ''
         self._tsmap        = gammalib.GSkyMap() # Empty sky map
-        self._statusmap  = gammalib.GSkyMap() # Empty sky map
+        self._statusmap    = gammalib.GSkyMap() # Empty sky map
         self._maps         = []
         self._mapnames     = []
         self._merged_files = []
@@ -148,7 +148,7 @@ class cstsmapmerge(ctools.cscript):
         for hdu in fits:
             # Leave out primary and status extension
-            if hdu.extname() != "IMAGE" and hdu.extname() != 'STATUS MAP':
+            if hdu.extname() != 'IMAGE' and hdu.extname() != 'STATUS MAP':
                 # Add present maps
                 skymap = gammalib.GSkyMap()
@@ -185,7 +185,7 @@ class cstsmapmerge(ctools.cscript):
         for hdu in fits:
             # Leave out primary and status extension
-            if hdu.extname() != "IMAGE" and hdu.extname() != 'STATUS MAP':
+            if hdu.extname() != 'IMAGE' and hdu.extname() != 'STATUS MAP':
                 # Add present maps
                 skymap = gammalib.GSkyMap()
@@ -197,9 +197,9 @@ class cstsmapmerge(ctools.cscript):
         # Compare size of maps
         if not len(add_maps) == len(self._maps):
-            msg = "Cannot merge map \""+fitsfile+"\" into map \""+\
-                  self._in_filename+"\" since the number of parameters "\
-                  "between both maps is different."
+            msg = 'Cannot merge map "'+fitsfile+'" into map "'+\
+                  self._in_filename+'" since the number of parameters ' +\
+                  'between both maps is different.'
             raise RuntimeError(msg) 
         # Loop over bins    
@@ -210,10 +210,10 @@ class cstsmapmerge(ctools.cscript):
                 # Raise exception if this bin has already been computed
                 if self._statusmap[i] > -0.5 and not self._overwrite:
-                    msg = "Attempt to merge bin which apparently has "+\
-                          "already been merged. File \""+fitsfile+"\" "+\
-                          "contains already merged bins. Set hidden "+\
-                          "parameter \"overwrite=yes\" to avoid this error."
+                    msg = 'Attempt to merge bin which apparently has '+\
+                          'already been merged. File "'+fitsfile+'" '+\
+                          'contains already merged bins. Set hidden '+\
+                          'parameter "overwrite=yes" to avoid this error.'
                     raise RuntimeError(msg)
                 # Copy TS values
@@ -240,7 +240,7 @@ class cstsmapmerge(ctools.cscript):
             Number of pixels for which TS has been computed.
         """
         # Get status map for this file
-        status = gammalib.GSkyMap(fitsfile+"[STATUS MAP]")
+        status = gammalib.GSkyMap(fitsfile+'[STATUS MAP]')
         # Count number of pixels with TS status set
         count = 0
@@ -265,10 +265,10 @@ class cstsmapmerge(ctools.cscript):
         self._get_parameters()
         # Write header into logger
-        self._log_header1(gammalib.TERSE, "Merge TS maps")
+        self._log_header1(gammalib.TERSE, 'Merge TS maps')
         # Initialise file to start with
-        file0              = ""
+        file0              = ''
         self._merged_files = []
         # Test files for the entry status map. Use the first one to appear
@@ -279,7 +279,7 @@ class cstsmapmerge(ctools.cscript):
             if not gammalib.GFilename(fitsfile).is_fits():
                 # Log message
-                self._log_value(gammalib.EXPLICIT, "Skip file", fitsfile +" (not a FITS file)")
+                self._log_value(gammalib.EXPLICIT, 'Skip file', fitsfile +' (not a FITS file)')
                 # Continue
                 continue
@@ -297,8 +297,8 @@ class cstsmapmerge(ctools.cscript):
                 file0 = fitsfile
                 # Log message
-                pix_info = " (%d TS pixels computed)" % self._get_number_of_ts_pixels(fitsfile)
-                self._log_value(gammalib.TERSE, "Initial TS map file", fitsfile + pix_info)
+                pix_info = ' (%d TS pixels computed)' % self._get_number_of_ts_pixels(fitsfile)
+                self._log_value(gammalib.TERSE, 'Initial TS map file', fitsfile + pix_info)
                 # Leave loop
                 break
@@ -310,16 +310,18 @@ class cstsmapmerge(ctools.cscript):
                 fits.close()  
                 # Info message             
-                self._log_value(gammalib.EXPLICIT, "Skip file", fitsfile +" (no \"STATUS MAP\" extension)")
+                self._log_value(gammalib.EXPLICIT, 'Skip file', fitsfile +' (no "STATUS MAP" extension)')
                 # Continue
                 continue
         # Signal if no suitable file was found
-        if file0 == "":
-            self._log_string(gammalib.TERSE, "None of the provided files seems to be a sliced "
-                          "TS map file (none has a \"STATUS MAP\" "
-                          "extension).\n")
+        if file0 == '':
+            msg = 'None of the provided files seems to be a sliced ' + \
+                          'TS map file (none has a "STATUS MAP" ' + \
+                          'extension)'
+            raise RuntimeError(msg)
+                          
         # ... otherwise merge files
         else:
@@ -341,7 +343,7 @@ class cstsmapmerge(ctools.cscript):
                 # Skip if file is not FITS
                 if not gammalib.GFilename(fitsfile).is_fits():
-                    self._log_value(gammalib.EXPLICIT, "Skip file", fitsfile +" (not a FITS file)")
+                    self._log_value(gammalib.EXPLICIT, 'Skip file', fitsfile +' (not a FITS file)')
                     continue   
                 # Open FITS file
@@ -350,15 +352,15 @@ class cstsmapmerge(ctools.cscript):
                 # Skip if file does not contain status map
                 if not fits.contains('STATUS MAP'):
                     fits.close()
-                    self._log_value(gammalib.EXPLICIT, "Skip file", fitsfile +" (no \"STATUS MAP\" extension)")
+                    self._log_value(gammalib.EXPLICIT, 'Skip file', fitsfile +' (no "STATUS MAP" extension)')
                     continue
                 # Close FITS file
                 fits.close()
                 # Logging
-                pix_info = " (%d TS pixels computed)" % self._get_number_of_ts_pixels(fitsfile)
-                self._log_value(gammalib.TERSE, "Merge TS map file", fitsfile + pix_info)
+                pix_info = ' (%d TS pixels computed)' % self._get_number_of_ts_pixels(fitsfile)
+                self._log_value(gammalib.TERSE, 'Merge TS map file', fitsfile + pix_info)
                 # Merge TS map
                 self._merge_ts_map(fitsfile)
@@ -375,7 +377,7 @@ class cstsmapmerge(ctools.cscript):
         """
         # Get output filename in case it was not read ahead
-        outmap = self["outmap"].filename()
+        outmap = self['outmap'].filename()
         # Create FITS file
         fits = gammalib.GFits()
@@ -402,29 +404,29 @@ class cstsmapmerge(ctools.cscript):
                 bins_merged += 1
         # Set boolean if map is fully done
-        done = (bins_merged == self._tsmaps.npix())
+        done = (bins_merged == self._tsmap.npix())
         # Log summary header
-        self._log_header2(gammalib.TERSE, "Merging Summary")
+        self._log_header2(gammalib.TERSE, 'Merging Summary')
-        # Write status map if we are not done yet
+        # Write computing status if we are not done yet
         if not done:
             # Log bins that were computed and those that were missing
-            self._log_value(gammalib.TERSE, "TS map bins", self._tsmap.npix())
-            self._log_value(gammalib.TERSE, "Bins computed and merged", bins_merged )
-            self._log_value(gammalib.TERSE, "Bins missing", self._tsmap.npix() - bins_merged)
+            self._log_value(gammalib.TERSE, 'TS map bins', self._tsmap.npix())
+            self._log_value(gammalib.TERSE, 'Bins computed and merged', bins_merged )
+            self._log_value(gammalib.TERSE, 'Bins missing', self._tsmap.npix() - bins_merged)
         else:
             # Write success message
-            self._log_string(gammalib.TERSE, "TS map was fully computed and successfully merged") 
+            self._log_string(gammalib.TERSE, 'TS map was fully computed and successfully merged') 
         # Write header
-        self._log_header1(gammalib.TERSE, "Save TS map")
+        self._log_header1(gammalib.TERSE, 'Save TS map')
         # Log filename
-        self._log_value(gammalib.TERSE, "TS map file", outmap.url())
+        self._log_value(gammalib.TERSE, 'TS map file', outmap.url())
         # Save FITS file
         fits.saveto(outmap, self._clobber()) 
@@ -433,7 +435,7 @@ class cstsmapmerge(ctools.cscript):
         if self._delete and done:
             # Log header
-            self._log_header2(gammalib.TERSE, "Delete slices TS map files")
+            self._log_header2(gammalib.TERSE, 'Delete slices TS map files')
             # Loop over sliced files
             for filename in self._merged_files:
@@ -442,7 +444,7 @@ class cstsmapmerge(ctools.cscript):
                 os.remove(filename)
                 # Log message about file
-                self._log_value(gammalib.TERSE,"Deleted input file", filename)
+                self._log_value(gammalib.TERSE,'Deleted input file', filename)
         # Return
         return
@@ -21,7 +21,10 @@ Gaussian sigma).
 cttsmap generates a FITS file comprising a sky map of TS values followed by 
 one extension per free parameter that contains sky maps of the fitted 
-parameter values.
+parameter values. To save computation time, the errors are not computed by
+default. Specifying the hidden parameter ``errors=yes`` will switch on the
+error computation. Accordingly, for each free parameter, the errors are also
+stored in a separate sky map.
 General parameters
@@ -56,6 +59,9 @@ General parameters
 ``outmap [file]``
     Output Test Statistic map file.
+    
+``(errors = no) [boolean]``
+    Compute and store parameter errors?
 ``(usepnt = no) [boolean]``
     Use CTA pointing direction for map centre instead of xref/yref parameters?
@@ -173,26 +173,20 @@ cttsmap&amp; cttsmap::operator=(const cttsmap&amp; app)
  ==========================================================================*/
 /***********************************************************************//**
- * @brief Clear cttsmap tool
- *
- * Clears cttsmap tool.
+ * @brief Clear instance
  ***************************************************************************/
 void cttsmap::clear(void)
 {
     // Free members
     free_members();
     this->ctool::free_members();
-
-    // Clear base class (needed to conserve tool name and version)
-    this->GApplication::clear();
+    this->GApplication::free_members();
     // Initialise members
+    this->GApplication::init_members();
     this->ctool::init_members();
     init_members();
-    // Write header into logger
-    log_header();
-
     // Return
     return;
 }
@@ -215,12 +209,9 @@ void cttsmap::run(void)
     // Get task parameters
     get_parameters();
-<<<<<<< HEAD
-=======
     // Write parameters into logger
     log_parameters(TERSE);
->>>>>>> b7dd76b35e6d2fd74ff010258ffa7606e072dbf5
     // Set energy dispersion flag for all CTA observations and save old
     // values in save_edisp vector
     std::vector<bool> save_edisp;
@@ -234,22 +225,16 @@ void cttsmap::run(void)
     }
     // Write observation(s) into logger
-    if (logTerse()) {
-        log << std::endl;
-        if (m_obs.size() > 1) {
-            log.header1("Observations");
-        }
-        else {
-            log.header1("Observation");
-        }
-        log << m_obs << std::endl;
+    if (m_obs.size() > 1) {
+    	log_header1(TERSE, "Observations");
+    }
+    else {
+    	log_header1(TERSE, "Observation");
     }
+    log_string(TERSE, m_obs.print());
     // Write header
-    if (logTerse()) {
-        log << std::endl;
-        log.header1("Initialise TS map");
-    }
+    log_header1(TERSE, "Initialise TS map");
     // Get bins to be computed
@@ -276,10 +261,7 @@ void cttsmap::run(void)
 	if (m_logL0 == 0.0) {
         // Write header
-        if (logTerse()) {
-            log << std::endl;
-            log.header1("Compute NULL Hypothesis for TS computation");
-        }
+		log_header1(TERSE, "Compute NULL Hypothesis for TS computation");
 		// Compute likelihood without the test source
 		m_obs.models(models);
@@ -289,10 +271,7 @@ void cttsmap::run(void)
 	}
     // Write header
-    if (logTerse()) {
-        log << std::endl;
-        log.header1("Generate TS map");
-    }
+    log_header1(TERSE, "Generate TS map");
     // Loop over grid positions
     for (int i = binmin; i < binmax; ++i) {
@@ -301,12 +280,9 @@ void cttsmap::run(void)
     	GSkyDir bincentre = m_tsmap.inx2dir(i);
     	// Header for verbose logging
-    	if (logExplicit()) {
-            std::string msg = "Computing TS for bin number "+gammalib::str(i)+
-                              " at "+bincentre.print();
-            log << std::endl;
-            log.header2(msg);
-    	}
+    	std::string msg = "Computing TS for bin number "+gammalib::str(i)+
+    	                              " at "+bincentre.print();
+    	log_header2(EXPLICIT, msg);
     	// Add test source at current bin position
         if (m_testsource != NULL) {
@@ -340,7 +316,6 @@ void cttsmap::run(void)
     		log << " TS value ....: ";
             log << ts << std::endl;
     	}
-
         else if (logTerse()) {
     		log << "TS for bin number ";
             log << i;
@@ -445,7 +420,7 @@ void cttsmap::save(void)
         // Write Fit status map
         m_statusmap.write(fits);
-        fits[fits.size()-1]->extname("FITSTATUS");
+        fits[fits.size()-1]->extname("STATUS MAP");
         // Save FITS file
         fits.saveto(m_outmap, clobber());
@@ -465,10 +440,7 @@ void cttsmap::save(void)
 void cttsmap::publish(const std::string& name)
 {
     // Write header
-    if (logTerse()) {
-        log << std::endl;
-        log.header1("Publish TS map");
-    }
+	log_header1(TERSE, "Publish TS map");
     // Set default name is user name is empty
     std::string user_name(name);
@@ -477,10 +449,7 @@ void cttsmap::publish(const std::string&amp; name)
     }
     // Log filename
-    if (logTerse()) {
-        log << gammalib::parformat("Publish TS map");
-        log << user_name << std::endl;
-    }
+    log_value(TERSE, "Publish TS map", user_name);
     // Publish TS map
     m_tsmap.publish(user_name);
@@ -659,9 +628,6 @@ void cttsmap::get_parameters(void)
         m_outmap = (*this)["outmap"].filename();
     }
-    // Write parameters into logger
-    log_parameters(TERSE);
-
     // Return
     return;
 }
@@ -716,7 +682,8 @@ void cttsmap::init_maps(const GSkyMap&amp; map)
 				// Store name of parameter error
 				m_mapnames.push_back("e_" + (*m_testsource)[i].name());
-            }
+
+            } // endif: error calculation was requested
         } // endfor: looped over all model parameters
@@ -27,7 +27,7 @@ bkgcube,   f, a, &quot;NONE&quot;,,, &quot;Input background cube file (only needed for stacked 
 caldb,     s, a, "prod2",,, "Calibration database"
 irf,       s, a, "South_0.5h",,, "Instrument response function"
 edisp,     b, h, no,,, "Apply energy dispersion?"
-errors,    b, h, no,,, "Compute/store parameter errors?"
+errors,    b, h, no,,, "Compute and store parameter errors?"
 outmap,    f, a, "tsmap.fits",,, "Output Test Statistic map file"
 #
@@ -21,6 +21,8 @@
 import gammalib
 import cscripts
 from testing import test
+from __builtin__ import False
+from Carbon.Aliases import false
 # ================================== #
@@ -98,7 +100,16 @@ class Test(test):
         # Check if execution failed
         self.test_assert(self._execute(cmd) != 0,
              'Check invalid input file when executed from command line')
-
+        
+        # Setup cstsmapmerge command
+        cmd = cstsmapmerge+' inmaps="'+self._events+' '+self._model+'"'+ \
+                           ' outmap="cstsmapmerge_cmd3.fits"'+ \
+                           ' logfile="cstsmapmerge_cmd3.log" chatter=4'
+                           
+                # Check if execution failed
+        self.test_assert(self._execute(cmd) != 0,
+             'Check input file that are no TS maps when executed from command line')
+              
         # Return
         return
@@ -160,10 +171,10 @@ class Test(test):
         # Check merged TS map
         self._check_result_file('cstsmapmerge_py4.fits')
-
-        # Set-up space-separated cstsmapmerge with invalid files
+        
+        # Set-up cstsmapmerge with incomplete file list
         merge = cscripts.cstsmapmerge()
-        merge['inmaps']  = self._events+' '+self._model
+        merge['inmaps']  = self._map0 + ' ' + self._events
         merge['outmap']  = 'cstsmapmerge_py5.fits'
         merge['logfile'] = 'cstsmapmerge_py5.log'
         merge['chatter'] = 4
@@ -172,13 +183,13 @@ class Test(test):
         merge.execute()
         # Check merged TS map
-        #self._check_result_file('cstsmapmerge_py5.fits')
+        self._check_result_file('cstsmapmerge_py5.fits', False)
         # Return
         return
     # Check result file
-    def _check_result_file(self, filename):
+    def _check_result_file(self, filename, test_complete = True):
         """
         Check result file
         """
@@ -189,6 +200,7 @@ class Test(test):
         ts        = fits['Primary']
         prefactor = fits['Prefactor']
         index     = fits['Index']
+        status    = fits['STATUS MAP']
         # Check dimensions
         self.test_value(ts.naxis(), 2, 'Check for 2 dimensions')
@@ -200,6 +212,23 @@ class Test(test):
         self.test_value(index.naxis(), 2, 'Check for 2 dimensions')
         self.test_value(index.naxes(0), 2, 'Check for 2 pixels in X')
         self.test_value(index.naxes(1), 1, 'Check for 1 pixel in Y')
+        self.test_value(status.naxis(), 2, 'Check for 2 dimensions')
+        self.test_value(status.naxes(0), 2, 'Check for 2 pixels in X')
+        self.test_value(status.naxes(1), 1, 'Check for 1 pixel in Y')
+        
+        # Initialise flag if map is complete
+        done = True
+        
+        # Loop over status sky map
+        for pix in status:
+            
+            # Check if pix is larger than threshold
+            if pix < -0.5:
+                done = False
+                break
+        
+        # Test for completeness of merged map
+        self.test_assert(done == test_complete, 'Test if map was merged completely')           
         # Return
         return