diff --git a/cscripts/cstsmapmerge.py b/cscripts/cstsmapmerge.py index e17c24f..9b12c8c 100755 --- a/cscripts/cstsmapmerge.py +++ b/cscripts/cstsmapmerge.py @@ -30,17 +30,17 @@ import ctools # ================== # class cstsmapmerge(ctools.cscript): """ - Merge Test Statistic maps + Merge Test Statistic maps. """ # Constructor def __init__(self, *argv): """ - Constructor + Constructor. """ # Set name self._name = 'cstsmapmerge' - self._version = '1.2.0' + self._version = '1.1.0' # Initialise class members self._files = None @@ -64,7 +64,7 @@ class cstsmapmerge(ctools.cscript): # Private methods def _get_parameters(self): """ - Get parameters from parfile + Get parameters from parfile. """ # Get input maps string inmaps = self['inmaps'].string() @@ -118,7 +118,7 @@ class cstsmapmerge(ctools.cscript): if self._read_ahead(): self['outmap'].filename() - # Write input parameters into logger + # Write input parameters into logger self._log_parameters(gammalib.TERSE) # Return @@ -126,12 +126,7 @@ class cstsmapmerge(ctools.cscript): def _init_ts_map(self, fitsfile): """ - Initialise Test Statistic map - - Parameters - ---------- - fitsfile : str - FITS file to be merged + Initialise Test Statistic map. """ # Set filename self._in_filename = fitsfile @@ -169,12 +164,10 @@ class cstsmapmerge(ctools.cscript): def _merge_ts_map(self, fitsfile): """ - Merge TS map from FITS file into output TS map + Merge TS map from FITS file into output TS map. - Parameters - ---------- - fitsfile : str - FITS file to be merged + Args: + fitsfile: FITS file to be merged. """ # Open FITS file fits = gammalib.GFits(fitsfile) @@ -204,23 +197,23 @@ class cstsmapmerge(ctools.cscript): # Compare size of maps if not len(add_maps) == len(self._maps): - msg = ('Cannot merge map "%s" into map "%s" since the number of ' - 'parameters between both maps is different.' % - (fitsfile, self._in_filename)) + 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 for i in range(self._tsmap.npix()): # Consider only bins that have been computed - if add_statusmap[i] > 0.5: + if add_statusmap[i] > -0.5: # 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 "%s" contains already merged ' - 'bins. Set hidden parameter "overwrite=yes" to ' - 'avoid this error.' % fitsfile) + 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.' raise RuntimeError(msg) # Copy TS values @@ -238,17 +231,13 @@ class cstsmapmerge(ctools.cscript): def _get_number_of_ts_pixels(self, fitsfile): """ - Return number of pixels with TS values + Return number of pixels with TS values. - Parameters - ---------- - fitsfile : str - FITS file to be merged + Args: + fitsfile: FITS file to be merged. - Returns - ------- - counts : int - Number of pixels for which TS has been computed + Returns: + Number of pixels for which TS has been computed. """ # Get status map for this file status = gammalib.GSkyMap(fitsfile+'[STATUS MAP]') @@ -256,7 +245,7 @@ class cstsmapmerge(ctools.cscript): # Count number of pixels with TS status set count = 0 for pix in status: - if pix > 0.5: + if pix > -0.5: count += 1 # Return number of pixels @@ -266,7 +255,7 @@ class cstsmapmerge(ctools.cscript): # Public methods def run(self): """ - Run the script + Run the script. """ # Switch screen logging on in debug mode if self._logDebug(): @@ -276,9 +265,7 @@ class cstsmapmerge(ctools.cscript): self._get_parameters() # Write header into logger - if self._logTerse(): - self._log('\n') - self._log.header1('Merge TS maps') + self._log_header1(gammalib.TERSE, 'Merge TS maps') # Initialise file to start with file0 = '' @@ -290,8 +277,11 @@ class cstsmapmerge(ctools.cscript): # Skip file if it's not a FITS file if not gammalib.GFilename(fitsfile).is_fits(): - self._log_value(gammalib.EXPLICIT, 'Skip file', - '%s (not a FITS file)' % fitsfile) + + # Log message + self._log_value(gammalib.EXPLICIT, 'Skip file', fitsfile +' (not a FITS file)') + + # Continue continue # Open FITS file @@ -299,29 +289,40 @@ class cstsmapmerge(ctools.cscript): # If file contains a status map then use it if fits.contains('STATUS MAP'): + + # Close FITS file fits.close() + + # Set filename to initial filename file0 = fitsfile - if self._logTerse(): - self._log(gammalib.parformat('Initial TS map file')) - self._log(fitsfile) - count = self._get_number_of_ts_pixels(fitsfile) - self._log(' (%d TS pixels computed)' % count) - self._log('\n') + + # 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) + + # Leave loop break # ... otherwise signal that file is useless else: - fits.close() - self._log_value(gammalib.EXPLICIT, 'Skip file', - '%s (no "STATUS MAP" extension)' % fitsfile) + + # Close fits file + fits.close() + + # Info message + self._log_value(gammalib.EXPLICIT, 'Skip file', fitsfile +' (no "STATUS MAP" extension)') + + # Continue continue # Signal if no suitable file was found if file0 == '': - if self._logTerse(): - self._log('None of the provided files seems to be a sliced ' - 'TS map file (none has a "STATUS MAP" extension).\n') - + 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: @@ -342,9 +343,8 @@ class cstsmapmerge(ctools.cscript): # Skip if file is not FITS if not gammalib.GFilename(fitsfile).is_fits(): - self._log_value(gammalib.EXPLICIT, 'Skip file', - '%s (not a FITS file)' % fitsfile) - continue + self._log_value(gammalib.EXPLICIT, 'Skip file', fitsfile +' (not a FITS file)') + continue # Open FITS file fits = gammalib.GFits(fitsfile) @@ -352,20 +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', - '%s (no "STATUS MAP" extension)' % fitsfile) + self._log_value(gammalib.EXPLICIT, 'Skip file', fitsfile +' (no "STATUS MAP" extension)') continue # Close FITS file fits.close() # Logging - if self._logTerse(): - self._log(gammalib.parformat('Merge TS map file')) - self._log(fitsfile) - count = self._get_number_of_ts_pixels(fitsfile) - self._log(' (%d TS pixels computed)' % count) - self._log('\n') + 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) @@ -378,16 +373,11 @@ class cstsmapmerge(ctools.cscript): def save(self): """ - Save TS map and remove slices if requested + Save TS map and remove slices if requested. """ - # Write header - self._log_header1(gammalib.TERSE, 'Save TS map') # Get output filename in case it was not read ahead outmap = self['outmap'].filename() - - # Log filename - self._log_value(gammalib.TERSE, 'TS map file', outmap.url()) # Create FITS file fits = gammalib.GFits() @@ -403,33 +393,65 @@ class cstsmapmerge(ctools.cscript): for i in range(len(self._mapnames)): fits[i+1].extname(self._mapnames[i]) + # Treat status map separately + self._statusmap.write(fits) + fits[fits.size()-1].extname('STATUS MAP') + # Check if map is fully done - done = True + bins_merged = 0 for pix in self._statusmap: - if pix < 0.5: - done = False - break + if pix > -0.5: + bins_merged += 1 + + # Set boolean if map is fully done + done = (bins_merged == self._tsmap.npix()) - # Write status map if we are not done yet + # Log summary header + self._log_header2(gammalib.TERSE, 'Merging Summary') + + # Write computing status if we are not done yet if not done: - self._statusmap.write(fits) - fits[fits.size()-1].extname('STATUS MAP') + + # 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) + + else: + + # Write success message + self._log_string(gammalib.TERSE, 'TS map was fully computed and successfully merged') + + # Write header + self._log_header1(gammalib.TERSE, 'Save TS map') + + # Log filename + self._log_value(gammalib.TERSE, 'TS map file', outmap.url()) # Save FITS file - fits.saveto(outmap, self._clobber()) + fits.saveto(outmap, self._clobber()) - # Delete TS input maps if requested - if self._delete: + # Delete TS input maps if requested and map is fully done + if self._delete and done: + + # Log header + self._log_header2(gammalib.TERSE, 'Delete slices TS map files') + + # Loop over sliced files for filename in self._merged_files: + + # Remove file os.remove(filename) - self._log_value(gammalib.TERSE, 'Deleted input file', filename) + + # Log message about file + self._log_value(gammalib.TERSE,'Deleted input file', filename) # Return return def execute(self): """ - Execute the script + Execute the script. """ # Open logfile self.logFileOpen() diff --git a/cscripts/cstsmapsplit.par b/cscripts/cstsmapsplit.par index 4736f66..a54e0ec 100644 --- a/cscripts/cstsmapsplit.par +++ b/cscripts/cstsmapsplit.par @@ -27,6 +27,7 @@ bkgcube, f, a, NONE,,, "Input background cube file (only needed for stacke 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?" outmap, f, a, tsmap.fits,,, "Output Test Statistic map file" bins_per_job, i, a, 5,,, "Number of TS map bins per task" compute_null, b, a, yes,,,"Compute null hypothesis first?" diff --git a/doc/source/users/reference_manual/cttsmap.rst b/doc/source/users/reference_manual/cttsmap.rst index aa26737..401125c 100644 --- a/doc/source/users/reference_manual/cttsmap.rst +++ b/doc/source/users/reference_manual/cttsmap.rst @@ -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? diff --git a/src/cttsmap/cttsmap.cpp b/src/cttsmap/cttsmap.cpp index ec8616d..7372bab 100644 --- a/src/cttsmap/cttsmap.cpp +++ b/src/cttsmap/cttsmap.cpp @@ -173,9 +173,7 @@ cttsmap& cttsmap::operator=(const cttsmap& app) ==========================================================================*/ /***********************************************************************//** - * @brief Clear cttsmap tool - * - * Clears cttsmap tool. + * @brief Clear instance ***************************************************************************/ void cttsmap::clear(void) { @@ -194,9 +192,6 @@ void cttsmap::clear(void) this->ctlikelihood::init_members(); init_members(); - // Write header into logger - log_header(); - // Return return; } @@ -290,6 +285,14 @@ void cttsmap::run(void) // Optimize observation container m_obs.optimize(m_opt); + // Compute errors if necessary + if (m_errors) { + m_obs.errors(m_opt); + } + + // Get status of optimization + int status = m_opt.status(); + // Retrieve the Likelihood value double logL1 = -(m_opt.value()); @@ -314,11 +317,27 @@ void cttsmap::run(void) // Extract fitted test source parameters for (int j = 0; j < m_mapnames.size(); ++j) { - m_maps[j](i) = (*testsource)[m_mapnames[j]].value(); - } - // Set status of bin to true - m_statusmap(i) = 1.0; + // If map name start with "e_" then set the parameter error ... + if (m_mapnames[j].substr(0,2) == "e_") { + + // Get parameter name by removing the error prefix + std::string parname = m_mapnames[j].substr(2, m_mapnames[j].size()); + + // Get parameter error + m_maps[j](i) = (*testsource)[parname].error(); + + } + + // ... otherwise set the parameter value + else { + m_maps[j](i) = (*testsource)[m_mapnames[j]].value(); + } + + } // endfor: looped over all maps + + // Set Fit status of the bin + m_statusmap(i) = status; // Remove model from container models.remove(m_srcname); @@ -371,11 +390,9 @@ void cttsmap::save(void) fits[i+1]->extname(m_mapnames[i]); } - // Add computation log if not all bins are computed - if (m_binmin != -1 || m_binmax != -1) { - m_statusmap.write(fits); - fits[m_mapnames.size()+1]->extname("STATUS MAP"); - } + // Write Fit status map + m_statusmap.write(fits); + fits[fits.size()-1]->extname("STATUS MAP"); // Save FITS file fits.saveto(m_outmap, clobber()); @@ -437,6 +454,7 @@ void cttsmap::init_members(void) m_outmap.clear(); m_apply_edisp = false; m_publish = false; + m_errors = false; // Initialise protected members m_binmin = -1; @@ -465,6 +483,7 @@ void cttsmap::copy_members(const cttsmap& app) m_outmap = app.m_outmap; m_apply_edisp = app.m_apply_edisp; m_publish = app.m_publish; + m_errors = app.m_errors; // Copy protected members m_binmin = app.m_binmin; @@ -516,7 +535,7 @@ void cttsmap::get_parameters(void) { // Setup observations from "inobs" parameter setup_observations(m_obs); - + // Setup models from "inmodel" parameter setup_models(m_obs, (*this)["srcname"].string()); @@ -545,6 +564,9 @@ void cttsmap::get_parameters(void) // Create sky map based on task parameters GSkyMap map = create_map(m_obs); + // Check if errors should be computed + m_errors = (*this)["errors"].boolean(); + // Initialise maps from user parameters init_maps(map); @@ -585,14 +607,17 @@ void cttsmap::init_maps(const GSkyMap& map) m_tsmap = GSkyMap(map); // Initialise map information - m_statusmap.clear(); + m_statusmap.clear(); + + // Create status map and progress map + m_statusmap = GSkyMap(map); - // Create status map - m_statusmap = GSkyMap(map); + // Set status map to invalid value -1 to signal that bin wasn't computed yet + m_statusmap = -1.0; // Loop over all model parameters for (int i = 0; i < m_testsource->size(); i++) { - + // Skip fixed parameters if ((*m_testsource)[i].is_fixed()) { continue; @@ -604,7 +629,18 @@ void cttsmap::init_maps(const GSkyMap& map) // Store parameter name m_mapnames.push_back((*m_testsource)[i].name()); - + + // Add parameter errors if requested + if (m_errors) { + + // Add sky map for fit error of the free parameter + m_maps.push_back(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 // Return diff --git a/src/cttsmap/cttsmap.hpp b/src/cttsmap/cttsmap.hpp index 7f2e12c..c4995b7 100644 --- a/src/cttsmap/cttsmap.hpp +++ b/src/cttsmap/cttsmap.hpp @@ -82,19 +82,20 @@ protected: void init_maps(const GSkyMap& map); // User parameters - std::string m_srcname; //!< Name of source which is moved around - GFilename m_outmap; //!< Output TS map file name - bool m_apply_edisp; //!< Apply energy dispersion? - bool m_publish; //!< Publish TS map? + std::string m_srcname; //!< Name of source which is moved around + GFilename m_outmap; //!< Output TS map file name + bool m_apply_edisp; //!< Apply energy dispersion? + bool m_publish; //!< Publish TS map? + bool m_errors; //!< Compute and store parameter errors? // Parameters to control speed and job splitting - int m_binmin; //!< Map bin number from which computation should start - int m_binmax; //!< Map bin number where map computation should end - double m_logL0; //!< Likelihood value of null hypothesis + int m_binmin; //!< Map bin number from which computation should start + int m_binmax; //!< Map bin number where map computation should end + double m_logL0; //!< Likelihood value of null hypothesis // Protected members GSkyMap m_tsmap; //!< TS map - GSkyMap m_statusmap; //!< Map of computed bins + GSkyMap m_statusmap; //!< Map of optimizer fit status std::vector m_mapnames; //!< Names of free parameters std::vector m_maps; //!< Sky maps for each free parameter GModel* m_testsource; //!< Pointer to test source for TS computation diff --git a/src/cttsmap/cttsmap.par b/src/cttsmap/cttsmap.par index 6c23137..70db7fc 100644 --- a/src/cttsmap/cttsmap.par +++ b/src/cttsmap/cttsmap.par @@ -5,11 +5,11 @@ # ------------------------------------------------------------------------- # # # # File: cttsmap.par # -# Version: 1.1.0 # +# Version: 1.2.0 # # # # Author: Michael Mayer # -# michael.mayer@desy.de # -# DESY # +# michael.mayer@physik.hu-berlin.de # +# HU Berlin # # # # Purpose: Parameter file for the cttsmap application # # # @@ -30,6 +30,11 @@ edisp, b, h, no,,, "Apply energy dispersion?" outmap, f, a, "tsmap.fits",,, "Output Test Statistic map file" # +# TS map computation parameters +#============================== +errors, b, h, no,,, "Compute and store parameter errors?" + +# # Spatial binning parameters #=========================== usepnt, b, h, no,,, "Use pointing instead of xref/yref parameters?" @@ -43,7 +48,7 @@ proj, s, a, "CAR", AIT|AZP|CAR|MER|MOL|STG|TAN,, "Projection method" # # Parameters for splitting and speed purpose -#=========================== +#=========================================== binmin, i, h, -1,,, "First bin to compute" binmax, i, h, -1,,, "Last bin to compute" logL0, r, h, 0.0,,, "LogLikelihood value of null hypothesis" diff --git a/test/data/tsmap_0.fits b/test/data/tsmap_0.fits index 57f530b..8aafb23 100644 --- a/test/data/tsmap_0.fits +++ b/test/data/tsmap_0.fits diff --git a/test/data/tsmap_1.fits b/test/data/tsmap_1.fits index e905c5b..01141ae 100644 --- a/test/data/tsmap_1.fits +++ b/test/data/tsmap_1.fits diff --git a/test/test_cstsmapmerge.py b/test/test_cstsmapmerge.py index c1dc352..cc80521 100755 --- a/test/test_cstsmapmerge.py +++ b/test/test_cstsmapmerge.py @@ -98,7 +98,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 +169,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 +181,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 +198,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 +210,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 diff --git a/test/test_cttsmap.py b/test/test_cttsmap.py index 0e2edea..864490f 100755 --- a/test/test_cttsmap.py +++ b/test/test_cttsmap.py @@ -97,6 +97,20 @@ class Test(test): # Check if execution failed self.test_assert(self._execute(cmd) != 0, 'Check invalid input file when executed from command line') + + # Setup cttsmap command + cmd = cttsmap+' inobs="'+self._events+'"'+ \ + ' inmodel="'+self._model+'" srcname="Crab"'+ \ + ' caldb="'+self._caldb+'" irf="'+self._irf+'"'+ \ + ' errors=yes'+ \ + ' outmap="cttsmap_cmd1.fits"'+ \ + ' nxpix=5 nypix=5 binsz=0.02'+ \ + ' coordsys="CEL" proj="CAR" xref=83.63 yref=22.01'+ \ + ' logfile="cttsmap_cmd1.log" chatter=1' + + # Check if execution was successful + self.test_assert(self._execute(cmd) == 0, + 'Check successful execution from command line') # Setup cttsmap --help option cmd = cttsmap+' --help'