Added cstsmapmerge and included documentation and unit tests

Michael Mayer
1 parent db573ea0
Showing 11 changed files with 630 additions and 13 deletions
cscripts/Makefile.am
cscripts/csmodelinfo.py
cscripts/csmodelmerge.py
cscripts/csobsinfo.py
cscripts/cstsmapmerge.py
doc/source/reference_manual/csmodelmerge.rst
doc/source/reference_manual/cstsmapmerge.rst
doc/source/reference_manual/index.rst
test/data/tsmap_0.fits
test/data/tsmap_1.fits
test/test_cscripts.sh.in
 ## Process this file with automake to produce Makefile.in
-# Scripts to be distributed in bindir
-binscripts = tsmerge.py
-
 # Scripts to be distributed in $(pyexecdir)/cscripts
 pythonscripts = __init__.py \
                 cspull.py cssens.py cstsdist.py cscaldb.py csresmap.py \
                 cshessobs.py csobsdef.py csspec.py cslightcrv.py csobsinfo.py \
-                csinfo.py csmodelinfo.py csmodelmerge.py \
+                csinfo.py csmodelinfo.py csmodelmerge.py cstsmapmerge.py \
                 obsutils.py
 # Scripts for which symbolic links will be created in bindir
 linkscripts = cspull cssens cstsdist cscaldb csresmap \
               cshessobs csobsdef csspec cslightcrv csobsinfo csmodelmerge \
-              csinfo csmodelinfo
+              csinfo csmodelinfo cstsmapmerge
 # Distribute scripts
 EXTRA_DIST         = $(binscripts) $(pythonscripts)
@@ -499,7 +499,7 @@ class csmodelinfo(ctools.cscript):
 # ======================== #
 if __name__ == '__main__':
     """
-    Generates caldb entry.
+    Generates model information output.
     """
     # Create instance of application
     app = csmodelinfo(sys.argv)
@@ -116,7 +116,7 @@ class csmodelmerge(ctools.cscript):
         # Throw exception if input models cannot be decoded
         else:
-            msg = "Parameter \"inmodels\" must contain either an @asci file, a comma-separated or space-separated list of files or a wildcard string"
+            msg = "Parameter \"inmodels\" must contain either an @ASCII file, a comma-separated or space-separated list of files or a wildcard string"
             raise gammalib.GException.invalid_argument("csmodelmerge::get_parameters", msg) 
         # Get output file
@@ -205,7 +205,7 @@ class csmodelmerge(ctools.cscript):
 # ======================== #
 if __name__ == '__main__':
     """
-    Generates caldb entry.
+    Merges model XML files
     """
     # Create instance of application
     app = csmodelmerge(sys.argv)
@@ -452,7 +452,7 @@ class csobsinfo(ctools.cscript):
 # ======================== #
 if __name__ == '__main__':
     """
-    Generates caldb entry.
+    Generates information about an observation container.
     """
     # Create instance of application
     app = csobsinfo(sys.argv)
+#!/usr/bin/env python
+# ==========================================================================
+# Merges sliced TS map files into one map
+#
+# Copyright (C) 2015 Michael Mayer
+#
+# This program is free software: you can redistribute it and/or modify
+# it under the terms of the GNU General Public License as published by
+# the Free Software Foundation, either version 3 of the License, or
+# (at your option) any later version.
+#
+# This program is distributed in the hope that it will be useful,
+# but WITHOUT ANY WARRANTY; without even the implied warranty of
+# MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
+# GNU General Public License for more details.
+#
+# You should have received a copy of the GNU General Public License
+# along with this program.  If not, see <http://www.gnu.org/licenses/>.
+#
+# ==========================================================================
+import gammalib
+import ctools
+import sys
+import glob
+import os
+
+
+# =============== #
+# cshessobs class #
+# =============== #
+class cstsmapmerge(ctools.cscript):
+    """
+    This class merges an arbitray number of model container into one output model XML file
+    """
+    def __init__(self, *argv):
+        """
+        Constructor.
+        """
+        # Set name
+        self.name    = "cstsmapmerge"
+        self.version = "0.1.0"
+
+        # Make sure that parfile exists
+        file = self.parfile()
+
+        # Initialise application
+        if len(argv) == 0:
+            ctools.cscript.__init__(self, self.name, self.version)
+        elif len(argv) == 1:
+            ctools.cscript.__init__(self, self.name, self.version, *argv)
+        else:
+            raise TypeError("Invalid number of arguments given.")
+
+        # Set logger properties
+        self.log_header()
+        self.log.date(True)
+
+        # Return
+        return
+
+    def __del__(self):
+        """
+        Destructor.
+        """
+        # Return
+        return
+
+    def parfile(self):
+        """
+        Check if parfile exists. If parfile does not exist then create a
+        default parfile. This kluge avoids shipping the cscript with a parfile.
+        """
+
+        # Set parfile name
+        parfile = self.name+".par"
+
+        try:
+            pars = gammalib.GApplicationPars(parfile)
+        except:
+            # Signal that parfile was not found
+            print("Parfile \""+parfile+"\" not found. Create default parfile.")
+
+            # Create default parfile
+            pars = gammalib.GApplicationPars()
+            pars.append(gammalib.GApplicationPar("inmaps","s","a","tsmap.fits","","","Input TS map FITS files"))
+            pars.append(gammalib.GApplicationPar("outmap","f","a","NONE","","","Output TS map FITS file"))
+            pars.append(gammalib.GApplicationPar("overwrite","b","h","yes","","","Overwrite previously filled bins?"))
+            pars.append(gammalib.GApplicationPar("delete","b","h","no","","","Delete merged files?"))
+            pars.append_standard()
+            pars.append(gammalib.GApplicationPar("logfile","f","h","cstsmapmerge.log","","","Log filename"))
+            pars.save(parfile)
+
+        # Return
+        return
+
+    def get_parameters(self):
+        """
+        Get parameters from parfile and setup the observation.
+        """
+        # Get input models
+        self.inmodels = self["inmaps"].string()
+        
+        # Interpet inmodels string
+        # Handle ascii files
+        if "@" == self.inmodels[0]:
+            self.files = open(self.inmodels.replace("@","")).read().splitlines()  
+            
+        # Handle wild card strings
+        elif "*" in self.inmodels:
+            self.files = glob.glob(self.inmodels)
+        
+        # Handle space-separated list
+        elif " " in self.inmodels:
+            self.files = self.inmodels.split(" ")
+        
+        # Handle comma-separated list
+        elif "," in self.inmodels:
+            self.files = self.inmodels.split(",")
+            
+        # Throw exception if input models cannot be decoded
+        else:
+            msg = "Parameter \"inmmaps\" must contain either an @ASCII file, a comma-separated or space-separated list of files or a wildcard string"
+            sys.exit(msg)
+                    
+        if len(self.files) <= 1:
+            msg = "Need at least two files to start merging, "+str(len(self.files))+" given"
+            sys.exit(msg)
+        
+        # Get output file
+        self.outmap = self["outmap"].filename()
+        
+        # Get other parameters
+        self.overwrite = self["overwrite"].boolean()
+        self.delete = self["delete"].boolean()
+        
+        # Return
+        return
+ 
+ 
+    def init_map(self, fitsfile):
+        
+        self.in_filename = fitsfile
+        
+        # open FITS file
+        
+        fits = gammalib.GFits(fitsfile)
+        self.tsmap = gammalib.GSkyMap()
+        self.tsmap.read(fits[0])
+        
+        self.statusmap = gammalib.GSkyMap()
+        self.statusmap.read(fits["STATUS MAP"])
+        
+        # Get other maps 
+        self.maps = []
+        self.mapnames = []
+        
+        # Loop over extensions
+        for hdu in fits:
+            
+            # Leave out primary and status extension
+            if hdu.extname() != "IMAGE" and hdu.extname() != "STATUS MAP":
+                
+                # Add present maps
+                skymap = gammalib.GSkyMap()
+                skymap.read(hdu)
+                self.maps.append(skymap)
+                self.mapnames.append(hdu.extname())
+                
+        # Close FITS file       
+        fits.close()
+        
+        # Return
+        return
+
+    # Add map instances to each other
+    # in this way the TS maps get merged  
+    def add(self, fitsfile):
+        
+        # open FITS file
+        fits = gammalib.GFits(fitsfile)
+        add_tsmap = gammalib.GSkyMap()
+        add_tsmap.read(fits[0])
+        
+        add_statusmap = gammalib.GSkyMap()
+        add_statusmap.read(fits["STATUS MAP"])
+        
+        # Get other maps 
+        add_maps = []
+        
+        # Loop over extensions
+        for hdu in fits:
+            
+            # Leave out primary and status extension
+            if hdu.extname() != "IMAGE" and hdu.extname() != "STATUS MAP":
+                
+                # Add present maps
+                skymap = gammalib.GSkyMap()
+                skymap.read(hdu)
+                add_maps.append(skymap)
+                
+        # Close FITS file
+        fits.close()
+        
+        # Compare size of other maps
+        if not len(add_maps) == len(self.maps):
+            self.error("Maps cannot be added: Maps from \""+self.in_filename+"\" do not match from \""+fitsfile+"\"")
+
+        # Loop over bins    
+        for i in range(self.tsmap.npix()):
+            if add_statusmap[i] > 0.5:
+                
+                # throw 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."
+                    msg += " Set hidden parameter overwrite=yes to avoid this error"
+                    raise gammalib.GException("cstsmapmerge::add", msg)
+                
+                # Copy TS values
+                self.tsmap[i] = add_tsmap[i]
+                
+                # Copy status 
+                self.statusmap[i] = add_statusmap[i]
+                
+                # Loop over maps and copy entries
+                for j in range(len(self.maps)):
+                    self.maps[j][i] = add_maps[j][i]
+        
+        # Return
+        return
+
+    def run(self):
+        """
+        Run the script.
+        """
+        # Switch screen logging on in debug mode
+        if self.logDebug():
+            self.log.cout(True)
+
+        # Get parameters
+        self.get_parameters()
+        
+        # Write input parameters into logger
+        if self.logTerse():
+            self.log_parameters()
+            self.log("\n")
+            
+        if self.logTerse():
+            self.log("\n")
+            self.log.header1("Merging sliced TS map files")
+            self.log("\n")
+        
+        # Initialise file to start with
+        file0 =""
+        
+        self.added_files = []
+        
+        # Test files for the entry status map
+        # use the first one to appear useful
+        for fitsfile in self.files:    
+            if not gammalib.is_fits(fitsfile):
+                self.log("Skipping file \""+fitsfile+" (not a FITS file)\n")
+            fits = gammalib.GFits(fitsfile)
+            
+            if fits.contains("STATUS MAP"):
+                self.log("Using file \""+fitsfile+"\" as initialising map\n")
+                file0 = fitsfile
+                break
+            else:
+                if self.logExplicit():
+                    self.log("File \""+fitsfile+"\" has not \"STATUS MAP\" extension\n")
+        
+        # Print an error message if no suitable file was provided
+        if file0 == "":
+            if self.logTerse():
+                self.log("None of the provided file seems to be a sliced ts map file (no \"STATUS MAP\" extension)\n") 
+        
+        # Otherwise merge files
+        else:
+            
+            # copy file list
+            workfiles = self.files
+            
+            # Remove entry which will be used to initalise the map
+            workfiles.remove(file0)
+            
+            # Initialise map from first file
+            self.init_map(file0) 
+            
+            # Append to added files 
+            self.added_files.append(file0)
+            
+            # Loop over files
+            for fitsfile in workfiles:
+                
+                # Skip if file is not FITS
+                if not gammalib.is_fits(fitsfile):
+                    if self.logTerse():
+                        self.log("Skipping file \""+fitsfile+" (not a FITS file)\n") 
+                    continue   
+                
+                # Skip if file does not contain status map
+                fits = gammalib.GFits(fitsfile)
+                if not fits.contains("STATUS MAP"):
+                    if self.logTerse():
+                        self.log("Skipping file \""+fitsfile+" (not a FITS file)\n")    
+                    continue
+                
+                # Logging
+                if self.logTerse():
+                    self.log("Adding TS map slice \""+fitsfile+"\"\n")  
+                    
+                # Add file 
+                self.add(fitsfile)
+                
+                # Append to added files
+                self.added_files.append(fitsfile)
+                
+        # Return
+        return
+                
+            
+    def save(self):
+        """ 
+        Save merged tsmap file and remove slices if requested
+        """
+        
+        # Logging
+        if self.logExplicit():
+            self.log("\n")
+            self.log.header1("Saving merged TS map")
+            self.log("\n") 
+        
+        # Create FITS file
+        fits = gammalib.GFits()
+        
+        # Write TS map into primary
+        self.tsmap.write(fits)
+
+        # Loop over maps and write them to fits
+        for i in range(len(self.maps)):
+            self.maps[i].write(fits)
+        
+        # Set map names as extensions
+        for i in range(len(self.mapnames)):   
+            fits[i+1].extname(self.mapnames[i])
+
+        # Check if map is fully done
+        done = True
+        for pix in self.statusmap:
+            if pix < 0.5:
+                done = False
+                break
+        
+        # Write Status map if we are not done yet
+        if not done:
+            self.statusmap.write(fits)
+            fits[fits.size()-1].extname("STATUS MAP")
+        
+        # Save FITS file
+        fits.saveto(self.outmap, self.clobber())
+        if self.logExplicit():
+            self.log("Saved output map to \""+self.outmap+"\"\n")
+        if self.delete:
+            if self.logTerse():
+                self.log("\n")
+                self.log.header1("Deleting slices")
+                self.log("\n")
+            for filename in self.added_files:
+                if self.logExplicit():
+                    self.log("Deleting file \""+filename+"\"\n")
+                os.remove(filename)
+        
+        # Return
+        return
+
+            
+    def execute(self):
+        """
+        Execute the script.
+        """
+        # Run the script
+        self.run()
+
+        # Save ds9 file if required
+        self.save()
+
+        # Return
+        return    
+        
+
+# ======================== #
+# Main routine entry point #
+# ======================== #
+if __name__ == '__main__':
+    """
+    Merges ts map slices.
+    """
+    # Create instance of application
+    app = cstsmapmerge(sys.argv)
+    
+    # Open logfile
+    app.logFileOpen()
+
+    # Execute application
+    app.execute()
+.. _csmodelmerge:
+
+csmodelmerge
+======
+
+Merges an arbitrary number of model XML files into one file.
+
+
+Synopsis
+--------
+
+This script can be used to merge several model XML file into one output file.
+The files can be passed via the ``inmodels`` parameter. There are 4 options to specify
+a list of files that should be merged:
+- a space-separated list of file names (e.g. model1.xml model2.xml)
+- a comma-separated list of file names (e.g. model1.xml,model2.xml)
+- a wildcard string (e.g. mymodel*.xml)
+- an ASCII file containing the file names, one file per line (e.g. @mymodels.txt)
+
+Note that the "@" needs to be specified in case of an ASCII file name. It is also important to 
+note that models with the same name cannot be appended. The tool will throw and exception if
+the provided file contain models with the same name.
+
+
+General parameters
+------------------
+
+``inmodels [string]``
+    Input model XML files
+
+``outmodel [file]``
+    Output model file
+    
+    
+Standard parameters
+-------------------
+
+``(chatter = 2) [integer]``
+    Verbosity of the executable:
+     chatter = 0: no information will be logged
+     
+     chatter = 1: only errors will be logged
+     
+     chatter = 2: errors and actions will be logged
+     
+     chatter = 3: report about the task execution
+     
+     chatter = 4: detailed report about the task execution
+ 	 	 
+``(clobber = yes) [boolean]``
+    Specifies whether an existing output counts cube should be overwritten.
+ 	 	 
+``(debug = no) [boolean]``
+    Enables debug mode. In debug mode the executable will dump any log file output to the console.
+ 	 	 
+``(mode = ql) [string]``
+    Mode of automatic parameters (default is "ql", i.e. "query and learn").
+
+``(logfile = cssens.log) [filename]``
+    Log filename.
+
+
+Related tools or scripts
+------------------------
+
+.. _cstsmapmerge:
+
+cstsmapmerge
+======
+
+Merges slices of splitted tsmap computations.
+
+
+Synopsis
+--------
+In case a TS map has been computed using the option to split the calculations into
+several jobs, this tool is required to merge the output into one final map again.
+Running :doc:`cttsmap` with the hidden parameters ``binmin`` and ``binmax`` specified 
+will result in several sliced TS maps. This script merges them into one FITS file.
+The input slices can be passed via the ``inmaps`` parameter. There are 4 options to specify
+a list of files that should be merged:
+- a space-separated list of file names (e.g. tsmap1.fits tsmap2.fits)
+- a comma-separated list of file names (e.g. tsmap1.fits,tsmap2.fits)
+- a wildcard string (e.g. tsmap*.xml)
+- an ASCII file containing the file names, one file per line (e.g. @mymaps.txt)
+
+Note that the "@" needs to be specified in case of an ASCII file name. 
+
+General parameters
+------------------
+
+``inmaps [string]``
+    Input TS map FITS files
+
+``outmap [file]``
+    Output TS map FITS file
+    
+``(overwrite=yes) [boolean]``
+    Overwrite previously filled bins?
+
+``(delete=no) [boolean]``
+    Delete merged files?
+    
+    
+Standard parameters
+-------------------
+
+``(chatter = 2) [integer]``
+    Verbosity of the executable:
+     chatter = 0: no information will be logged
+     
+     chatter = 1: only errors will be logged
+     
+     chatter = 2: errors and actions will be logged
+     
+     chatter = 3: report about the task execution
+     
+     chatter = 4: detailed report about the task execution
+ 	 	 
+``(clobber = yes) [boolean]``
+    Specifies whether an existing output counts cube should be overwritten.
+ 	 	 
+``(debug = no) [boolean]``
+    Enables debug mode. In debug mode the executable will dump any log file output to the console.
+ 	 	 
+``(mode = ql) [string]``
+    Mode of automatic parameters (default is "ql", i.e. "query and learn").
+
+``(logfile = cssens.log) [filename]``
+    Log filename.
+
+
+Related tools or scripts
+------------------------
+
@@ -47,3 +47,5 @@ cscripts
    cstsdist --- Generates TS distribution <cstsdist>
    csobsinfo --- Shows detailled information about an observation container
    csmodelinfo --- Shows detailled information about a model container
+   csmodelmerge --- Merges several model containers into one file
+   cstsmapmerge --- Merges slices from ts map computations
@@ -117,7 +117,7 @@ $ECHO &quot; ok&quot;
 #
 # Test csmodelmerge (space-separated)
 # ==============
-$ECHO -n "Test csmodelmerge: "
+$ECHO -n "Test csmodelmerge (space-separated): "
 chmod +x @abs_top_srcdir@/cscripts/csmodelmerge.py
 @abs_top_srcdir@/cscripts/csmodelmerge.py \
   inmodels="data/crab.xml data/source1.xml" \
@@ -136,7 +136,7 @@ $ECHO &quot; ok&quot;
 #
 # Test csmodelmerge (comma-separated)
 # ==============
-$ECHO -n "Test csmodelmerge: "
+$ECHO -n "Test csmodelmerge (comma-separated): "
 chmod +x @abs_top_srcdir@/cscripts/csmodelmerge.py
 @abs_top_srcdir@/cscripts/csmodelmerge.py \
   inmodels="data/crab.xml,data/source1.xml" \
@@ -155,7 +155,7 @@ $ECHO &quot; ok&quot;
 #
 # Test csmodelmerge (wildcard string)
 # ==============
-$ECHO -n "Test csmodelmerge: "
+$ECHO -n "Test csmodelmerge (wildcard string): "
 chmod +x @abs_top_srcdir@/cscripts/csmodelmerge.py
 @abs_top_srcdir@/cscripts/csmodelmerge.py \
   inmodels="data/model_cube_background*.xml" \
@@ -175,7 +175,7 @@ $ECHO &quot; ok&quot;
 # Test csmodelmerge (ascii-file)
 # ==============
 ls data/model_cube_background* > models.txt
-$ECHO -n "Test csmodelmerge: "
+$ECHO -n "Test csmodelmerge (ascii-file): "
 chmod +x @abs_top_srcdir@/cscripts/csmodelmerge.py
 @abs_top_srcdir@/cscripts/csmodelmerge.py \
   inmodels="@models.txt" \
@@ -193,6 +193,84 @@ rm models.txt
 #
+# Test cstsmapmerge (comma-separated)
+# ==============
+$ECHO -n "Test cstsmapmerge (comma-separated): "
+chmod +x @abs_top_srcdir@/cscripts/cstsmapmerge.py
+@abs_top_srcdir@/cscripts/cstsmapmerge.py \
+  inmaps="data/tsmap_0.fits,data/tsmap_1.fits" \
+  outmap="merged_tsmap1.fits" > cstsmapmerge.log
+$ECHO -n "."
+if [ -s "merged_tsmap1.fits" ]
+then
+  $ECHO -n "."
+else
+  $ECHO " merged_tsmap1.fits file is not valid"
+  exit 1
+fi
+$ECHO " ok"
+
+
+#
+# Test cstsmapmerge (space-separated)
+# ==============
+$ECHO -n "Test cstsmapmerge (space-separated): "
+chmod +x @abs_top_srcdir@/cscripts/cstsmapmerge.py
+@abs_top_srcdir@/cscripts/cstsmapmerge.py \
+  inmaps="data/tsmap_0.fits data/tsmap_1.fits" \
+  outmap="merged_tsmap2.fits" > cstsmapmerge.log
+$ECHO -n "."
+if [ -s "merged_tsmap2.fits" ]
+then
+  $ECHO -n "."
+else
+  $ECHO " merged_tsmap2.fits file is not valid"
+  exit 1
+fi
+$ECHO " ok"
+
+
+#
+# Test cstsmapmerge (wildcard string)
+# ==============
+$ECHO -n "Test cstsmapmerge (wildcard string): "
+chmod +x @abs_top_srcdir@/cscripts/cstsmapmerge.py
+@abs_top_srcdir@/cscripts/cstsmapmerge.py \
+  inmaps="data/tsmap_*.fits" \
+  outmap="merged_tsmap3.fits" > cstsmapmerge.log
+$ECHO -n "."
+if [ -s "merged_tsmap3.fits" ]
+then
+  $ECHO -n "."
+else
+  $ECHO " merged_tsmap3.fits file is not valid"
+  exit 1
+fi
+$ECHO " ok"
+
+
+#
+# Test cstsmapmerge (ascii-file)
+# ==============
+ls data/tsmap_*.fits > tsmaps.txt
+$ECHO -n "Test cstsmapmerge  (ascii-file): "
+chmod +x @abs_top_srcdir@/cscripts/cstsmapmerge.py
+@abs_top_srcdir@/cscripts/cstsmapmerge.py \
+  inmaps="@tsmaps.txt" \
+  outmap="merged_tsmap4.fits" > cstsmapmerge.log
+$ECHO -n "."
+if [ -s "merged_tsmap4.fits" ]
+then
+  $ECHO -n "."
+else
+  $ECHO " merged_tsmap4.fits file is not valid"
+  exit 1
+fi
+$ECHO " ok"
+rm tsmaps.txt
+
+
+#
 # Test cspull
 # ===========
 $ECHO -n "Test cspull: "