Merge branch '1797-cstssplit' into integration

Jürgen Knödlseder
2 parents d8fcfdf3 b6c89800
Showing 8 changed files with 687 additions and 2 deletions
cscripts/Makefile.am
cscripts/__init__.py.in
cscripts/cstsmapsplit.par
cscripts/cstsmapsplit.py
doc/source/users/reference_manual/cstsmapsplit.rst
doc/source/users/reference_manual/index.rst
test/test_cstsmapsplit.py
test/test_python_cscripts.py
@@ -16,6 +16,7 @@ pythonscripts = __init__.py \
                 cssens.py \
                 csspec.py \
                 cstsdist.py \
+                cstsmapsplit.py \
                 cstsmapmerge.py \
                 csworkflow.py \
                 csfindobs.py \
@@ -38,6 +39,7 @@ parfiles = cscaldb.par \
            cssens.par \
            csspec.par \
            cstsdist.par \
+           cstsmapsplit.par \
            cstsmapmerge.par \
            csworkflow.par \
            csfindobs.par \
@@ -60,6 +62,7 @@ linkscripts = cscaldb \
               cssens \
               csspec \
               cstsdist \
+              cstsmapsplit \
               cstsmapmerge \
               csworkflow \
               csfindobs \
@@ -84,6 +87,7 @@ testscripts = tests/__init__.py \
               ../test/test_cssens.py \
               ../test/test_csspec.py \
               ../test/test_cstsdist.py \
+              ../test/test_cstsmapsplit.py \
               ../test/test_cstsmapmerge.py \
               ../test/test_csworkflow.py \
               ../test/test_csfindobs.py \
@@ -38,6 +38,7 @@ __all__ = [
     "cssens",
     "csspec",
     "cstsdist",
+    "cstsmapsplit",
     "cstsmapmerge",
     "csworkflow",
     "obsutils"
@@ -63,6 +64,7 @@ from cscripts.csroot2caldb import csroot2caldb
 from cscripts.cssens       import cssens
 from cscripts.csspec       import csspec
 from cscripts.cstsdist     import cstsdist
+from cscripts.cstsmapsplit import cstsmapsplit
 from cscripts.cstsmapmerge import cstsmapmerge
 from cscripts.csworkflow   import csworkflow
 from cscripts              import obsutils
+#############################################################################
+#                                                                           #
+#                       CTA Science Analysis Tools                          #
+#                                                                           #
+# ------------------------------------------------------------------------- #
+#                                                                           #
+#  File:      cstsmapsplit.par                                              #
+#  Version:   1.1.0                                                         #
+#                                                                           #
+#  Author:    Michael Mayer                                                 #
+#             michael.mayer@physik.hu-berlin.de                             #
+#             HU Berlin                                                     #
+#                                                                           #
+#  Purpose:   Parameter file for the cstsmapsplit application               #
+#                                                                           #
+#############################################################################
+#
+# 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"
+srcname,       s, a, Crab,,, "Test source name"
+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?"
+outmap,        f, a, tsmap.fits,,, "Output Test Statistic map file"
+bins_per_job,  i, a, 5,,, "Number of TS map bins per running job"
+compute_null,  b, a, yes,,,"Compute null hypothesis first?"
+outfile,       f, a, commands.txt,,, "ASCII file containing all commands"
+
+#
+# Spatial binning parameters
+#===========================
+usepnt,    b, h, no,,, "Use pointing instead of xref/yref parameters?"
+nxpix,     i, a, 200,,, "Size of the X axis in pixels"
+nypix,     i, a, 200,,, "Size of the Y axis in pixels"
+binsz,     r, a, 0.02,,, "Image scale (in degrees/pixel)"
+coordsys,  s, a, CEL, CEL|GAL,,"Coordinate system (CEL - celestial, GAL - galactic)"
+xref,      r, a, 83.63,0,360, "First coordinate of image center in degrees (RA or galactic l)"
+yref,      r, a, 22.01,-90,90, "Second coordinate of image center in degrees (DEC or galactic b)"
+proj,      s, a, CAR, AIT|AZP|CAR|MER|MOL|STG|TAN,, "Projection method"
+
+#
+# Standard parameters
+#====================
+publish, b, h, no,,, "Publish Test Statistic map on VO Hub?"
+chatter, i, h, 2,0,4, "Chattiness of output"
+clobber, b, h, yes,,, "Overwrite existing output files with new output files?"
+debug,   b, h, no,,, "Debugging mode activated"
+mode,    s, h, "ql",,, "Mode of automatic parameters"
+logfile, f, h, cstsmapsplit.log,,, "Log filename"
+#!/usr/bin/env python
+# ==========================================================================
+# Create commands to split TS map computation
+#
+# Copyright (C) 2016 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 sys
+import os
+import math
+import gammalib
+import ctools
+
+
+# ================== #
+# cstsmapsplit class #
+# ================== #
+class cstsmapsplit(ctools.cscript):
+    """
+    Generates commands to split TS map computation
+    
+    This class implements the creation of file containing commands to split
+    the computation of a TS map into several machines. It derives from
+    the ctools.cscript class which provides support for parameter files,
+    command line arguments, and logging. In that way the Python script
+    behaves just as a regular ctool. 
+    """
+    
+    # Consructor
+    def __init__(self, *argv):
+        """
+        Constructor
+        """
+
+        # Set name and version
+        self._name    = 'cstsmapsplit'
+        self._version = '1.1.0'
+
+        # Set data members
+        self._outmap       = gammalib.GFilename()
+        self._bins_per_job = 0
+        self._compute_null = False
+        self._outfile      = gammalib.GFilename()
+        self._map          = gammalib.GSkyMap()
+        self._cmd          = []
+        
+        # Initialise observation container from constructor arguments.
+        self._obs, argv = self._set_input_obs(argv)
+
+        # Initialise application by calling the appropriate class
+        # constructor.
+        self._init_cscript(argv)
+
+        # Return
+        return
+    
+    
+    # Private methods
+    def _get_parameters(self):
+        """
+        Get parameters from parfile and setup the observation
+        """
+         
+        # If there are no observations in container then get some ...
+        if self._obs.size() == 0:
+            self._obs = self._get_observations()
+            
+        # Set models if we have none
+        if self._obs.models().size() == 0:
+            self._obs.models(self['inmodel'].filename())
+        
+        # Get TS map paramters
+        self._map      = self._create_map(self._obs)        
+        self._srcname  = self['srcname'].string()      
+        self._outmap   = self['outmap'].filename()
+        
+        # Get additional parameters
+        self._bins_per_job = self['bins_per_job'].integer()
+        self._compute_null = self['compute_null'].boolean()
+        self._outfile      = self['outfile'].filename()
+        
+        # Write input parameters into logger
+        if self._logTerse():
+            self._log_parameters()
+            self._log('\n')
+        
+        # Return
+        return
+ 
+    
+    def _compute_null_hypothesis(self):
+        """
+        Compute null hypothesis
+        
+        Returns
+        -------
+        logl : float
+            Log-likelihood of null hypothesis
+        """
+        
+        # Store original models  
+        models_orig = self._obs.models()
+        
+        # Get models instance
+        models0 = self._obs.models() 
+        
+        # Remove test source   
+        models0.remove(self._srcname)
+        
+        # Initialise optimizer
+        opt = gammalib.GOptimizerLM()        
+        
+        # Set models for null hypothesis 
+        self._obs.models(models0)
+        
+        # Optimise observation container
+        self._obs.optimize(opt)
+        
+        # Set original models to container
+        self._obs.models(models_orig)
+              
+        # Return optimised log-likelihood of null hypothesis
+        return -(opt.value())
+    
+    
+    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 information into logger
+        if self._logTerse():
+            self._log('\n')
+            self._log.header1('Test source')
+            model = self._obs.models()[self._srcname]
+            self._log(str(model))
+            self._log('\n')
+        
+        # Set log-likelihood to zero
+        logL0 = 0.0
+        
+        # Pre-compute null hypothesis if requested
+        if self._compute_null:
+            
+            # Write information into logger
+            if self._logTerse():
+                self._log('\n')
+                self._log.header1('Compute null hypothesis')
+            
+            # Compute null hypothesis
+            logL0 = self._compute_null_hypothesis()
+            
+            # Write likelihood into logger
+            if self._logTerse():
+                self._log(gammalib.parformat('Source removed'))
+                self._log(self._srcname)
+                self._log('\n')
+                self._log(gammalib.parformat('Log-likelihood'))
+                self._log(repr(logL0))
+                self._log('\n')
+        
+        # Get parameters of TS map ctool
+        pars = gammalib.GApplicationPars('cttsmap.par')
+        
+        # Compute total number of jobs
+        nbins = self._map.npix()
+        njobs = int(math.ceil(float(nbins) / float(self._bins_per_job)) + 0.1)   
+        
+        # Set parameters to be skipped now, we will deal with them later
+        skip_pars = ['binmin', 'binmax', 'logL0', 'outmap', 'logfile']
+        
+        # Set tool name for computation
+        base_command = 'cttsmap'
+        
+        # Write information into logger
+        if self._logTerse():
+            self._log('\n')
+            self._log.header1('Create commands')
+            self._log(gammalib.parformat('Number of cttsmap calls'))
+            self._log(str(njobs))
+            self._log('\n')
+        
+        # Loop over TS map parameters
+        for par in pars:
+            
+            # Skip if we deal with them later
+            if par.name() in skip_pars:
+                continue
+            
+            # Skip if they need to be queried
+            # This way we ensure we pass only parameters
+            # that were queried before
+            if self[par.name()].is_query():
+                continue
+
+            # Set TS map parameter according to input from
+            # this script
+            par.value(self[par.name()].value())
+            
+            # Append command to set parameter
+            base_command += ' ' + par.name() + '=' + par.value()
+            
+        # Append null hypothesis parameter
+        base_command += ' logL0=' + repr(logL0)  
+        
+        # Set binning to start from zero
+        binmin = 0
+        binmax = 0
+        
+        # Clear command sequence
+        self._cmd = []
+        
+        # Loop over jobs and create commands
+        for job in range(njobs):
+            
+            # Set specific outmap file name
+            outmap = self._outmap.url().replace('.fits', '_'+str(job)+'.fits')
+            
+            # Set bin numbers to be computed in this job
+            binmin = binmax
+            binmax = binmin + self._bins_per_job
+            if binmax > nbins:
+                binmax = nbins           
+                
+            # Setup sliced command
+            sliced_command = '%s binmin=%s binmax=%s outmap=%s logfile=%s' % \
+                             (base_command, str(binmin), str(binmax),
+                              outmap, outmap.replace('.fits','.log'))
+            
+            # Append command to list of commands
+            self._cmd.append(sliced_command)
+        
+        # Write information into logger
+        if self._logExplicit():
+            self._log('\n')
+            self._log.header2('Commands')
+            for cmd in self._cmd:
+                self._log(cmd)
+                self._log('\n')
+            self._log('\n')
+
+        # Return
+        return
+    
+    
+    def execute(self):
+        """
+        Execute the script
+        """
+        
+        # Run the script
+        self.run()
+        
+        # Save results
+        self.save()
+        
+        # Return
+        return
+
+    
+    def save(self):
+        """
+        Save commands to ASCII file
+        """
+        
+        # Open file
+        f = open(self._outfile.url(), 'w')
+        
+        # Write commands to file
+        for cmd in self._cmd:
+            f.write(cmd + '\n')
+            
+        # Close file
+        f.close()
+        
+        # Return
+        return
+
+
+# ======================== #
+# Main routine entry point #
+# ======================== #
+if __name__ == '__main__':
+    
+    # Create instance of application
+    app = cstsmapsplit(sys.argv)
+    
+    # Execute application
+    app.execute()
+.. _cstsmapsplit:
+
+cstsmapsplit
+============
+
+Generate commands to split the Test Statistic map computation
+
+
+Synopsis
+--------
+
+This script generates a sequence of commands to split the Test Statistic (TS)
+map computation into a distinct number of :ref:`cttsmap` computations. The
+sequence of commands is written into an ASCII file, where each line in the
+file corresponds to the command for an individual :ref:`cttsmap` run. Each TS
+map computed in that way will be suffixed with the command number, and the
+map can be combined into a single TS map using the :ref:`cstsmapmerge` script.
+
+
+General parameters
+------------------
+
+``inobs [file]``
+    Input event list, counts cube or observation definition XML file.
+
+``inmodel [file]``
+    Input model XML file.
+
+``srcname [string]``
+    Name of source model for which the TS map should be computed.
+
+``expcube [file]``
+    Input exposure cube file (only needed for stacked analysis).
+
+``psfcube [file]``
+    Input PSF cube file (only needed for stacked analysis).
+
+``bkgcube [file]``
+    Input background cube file (only needed for stacked analysis).
+
+``caldb [string]``
+    Calibration database.
+
+``irf [string]``
+    Instrumental response function.
+
+``(edisp = no) [boolean]``
+    Applies energy dispersion to response computation.
+
+``outmap [file]``
+    Output Test Statistic map file.
+
+``bins_per_job [integer]``
+    Number of bins of the TS map that should be computed within one command
+
+``compute_null [boolean]``
+    Pre-computes the null hypothesis globally and passes the value to the
+    individual commands.  If the parameter is set to ``no``, the null
+    hypothesis will be computed in each run of :ref:`cttsmap` individually.
+    
+``outfile [file]``
+	Output ASCII file name where the commands are written to.
+
+``(usepnt = no) [boolean]``
+    Use CTA pointing direction for map centre instead of xref/yref parameters?
+ 	 	 
+``nxpix [integer]``
+    Number of map pixels in Right Ascension or Galactic longitude.
+ 	 	 
+``nypix [integer]``
+    Number of map pixels in Declination or Galactic latitude.
+ 	 	 
+``binsz [real]``
+    Map pixel size (in degrees/pixel).
+ 	 	 
+``coordsys <CEL|GAL> [string]``
+    Coordinate system (CEL - celestial, GAL - galactic).
+ 	 	 
+``xref [real]``
+    Right Ascension / Galactic longitude of map centre (J2000, in degrees).
+ 	 	 
+``yref [real]``
+    Declination / Galactic latitude of map centre (J2000, in degrees).
+ 	 	 
+``proj <AIT|AZP|CAR|MER|MOL|STG|TAN> [string]``
+    Projection method.
+
+
+
+Standard parameters
+-------------------
+
+``(publish = no) [boolean]``
+    Specifies whether the TS map should be published on VO Hub.
+
+``(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 TS map file 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 = cstsmapsplit.log) [string]``
+    Name of log file.
+
+
+Related tools or scripts
+------------------------
+
+:ref:`cttsmap`
+:ref:`cstsmapmerge`
+
@@ -49,8 +49,10 @@ cscripts
    csresmap --- Generates residual map <csresmap>
    cssens --- Computes CTA sensitivity <cssens>
    csspec --- Computes spectral points <csspec>
-   cstsdist --- Generates TS distribution <cstsdist>
-   cstsmapmerge --- Merges slices from ts map computations <cstsmapmerge>
+   cstsdist --- Generates Test Statistic distribution <cstsdist>
+   cstsmapsplit --- Creates commands to split the Test Statistic map computations <cstsmapsplit>
+   cstsmapmerge --- Merges slices from Test Statistic map computations <cstsmapmerge>
+   
  
  
 Scripts to manage an IACT database
+#! /usr/bin/env python
+# ==========================================================================
+# This scripts performs unit tests for the cstsmapsplit script.
+#
+# Copyright (C) 2016 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 os
+import gammalib
+import cscripts
+
+
+# =============================== #
+# Test class for cstssplit script #
+# =============================== #
+class Test(gammalib.GPythonTestSuite):
+    """
+    Test class for cstssplit script
+
+    This test class makes unit tests for the csstsplit script by using it
+    from the command line and from Python.
+    """
+    
+    # Constructor
+    def __init__(self):
+        """
+        Constructor
+        """
+        # Call base class constructor
+        gammalib.GPythonTestSuite.__init__(self)
+
+        # Return
+        return
+
+    # Set test functions
+    def set(self):
+        """
+        Set all test functions
+        """
+        # Set test name
+        self.name('cstsmapsplit')
+
+        # Append tests
+        self.append(self._test_cmd, 'Test cstsmapsplit on command line')
+        #self.append(self._test_python, 'Test cstsmapsplit from Python')
+
+        # Return
+        return
+
+    # Test cstsmapsplit on command line
+    def _test_cmd(self):
+        """
+        Test cstsmapsplit on the command line.
+        """
+        # Kluge to set the command (installed version has no README file)
+        if os.path.isfile("README.md"):
+            cstsmapsplit = "../cscripts/cstsmapsplit.py"
+        else:
+            cstsmapsplit = "cstsmapsplit"
+
+        # Setup cstsmapsplit command
+        cmd = cstsmapsplit+' inobs="selected_events.fits"'+ \
+                           ' inmodel="data/crab.xml"'+ \
+                           ' srcname="Crab"'+ \
+                           ' caldb="prod2"'+ \
+                           ' irf="South_0.5h"'+ \
+                           ' outmap="tsmap.fits"'+ \
+                           ' nxpix=5'+ \
+                           ' nypix=5'+ \
+                           ' binsz=0.05'+ \
+                           ' coordsys="CEL"'+ \
+                           ' xref=83.63'+ \
+                           ' yref=22.01'+ \
+                           ' proj="CAR"'+ \
+                           ' compute_null=yes'+ \
+                           ' bins_per_job=5'+ \
+                           ' outfile="tsmap_cmds1.txt"'
+
+
+        # Execute cstssplit, make sure we catch any exception
+        try:
+            rc = os.system(cmd+" >/dev/null 2>&1")
+        except:
+            pass
+
+        # Check if execution was successful
+        self.test_assert(rc == 0,
+                         "Successful cstsssplit execution on command line")
+
+        # Check observation definition XML file
+        self._check_cmdfile('tsmap_cmds1.txt')
+
+        # Setup cstsmapsplit command
+        cmd = cstsmapsplit+' inobs="selected_events.fits"'+ \
+                           ' inmodel="data/crab.xml"'+ \
+                           ' srcname="Does_not_exist"'+ \
+                           ' caldb="prod2"'+ \
+                           ' irf="South_0.5h"'+ \
+                           ' outmap="tsmap.fits"'+ \
+                           ' nxpix=5'+ \
+                           ' nypix=5'+ \
+                           ' binsz=0.05'+ \
+                           ' coordsys="CEL"'+ \
+                           ' xref=83.63'+ \
+                           ' yref=22.01'+ \
+                           ' proj="CAR"'+ \
+                           ' compute_null=yes'+ \
+                           ' bins_per_job=5'+ \
+                           ' outfile="tsmap_cmds2.txt"'
+                        
+
+        # Execute cstssplit, make sure we catch any exception
+        try:
+            rc = os.system(cmd+" >/dev/null 2>&1")
+        except:
+            pass
+
+        # Check if execution failed
+        self.test_assert(rc != 0,
+                         "Failure of cstssplit execution on command line")
+
+        # Return
+        return
+
+    # Test cstsmapsplit from Python
+    def _test_python(self):
+        """
+        Test cstsmapsplit from Python.
+        """
+        # Set-up cstsmapsplit
+        tsmapsplit = cscripts.cstsmapsplit()
+        tsmapsplit['inobs']        = 'obs_selected.xml'
+        tsmapsplit['inmodel']      = 'data/crab.xml'
+        tsmapsplit['xref']         = 83.6331
+        tsmapsplit['yref']         = 22.01
+        tsmapsplit['coordsys']     = 'CEL'
+        tsmapsplit['proj']         = 'CAR'
+        tsmapsplit['binsz']        = 0.05
+        tsmapsplit['nxpix']        = 5
+        tsmapsplit['nypix']        = 5
+        tsmapsplit['compute_null'] = False
+        tsmapsplit['bins_per_job'] = 5
+        tsmapsplit['outfile']      = 'tsmap_cmd_py1.txt'
+
+        # Run cstssplit script and save run list
+        tsmapsplit.logFileOpen()   # Make sure we get a log file
+        tsmapsplit.run()
+        tsmapsplit.save()
+
+        # Check observation definition XML file
+        self._check_cmdfile('tsmap_cmd_py1.txt')
+        
+        # Return
+        return
+
+    # Check observation definition XML file
+    def _check_cmdfile(self, filename):
+        """
+        Check ascii files containing commands
+        """
+        
+        # Create file name instance
+        fname = gammalib.GFilename(filename)
+        
+        # Check if execution was successful
+        self.test_assert(fname.exists(), 'Check of output filename exists')
+        
+        # Return
+        return
@@ -37,6 +37,7 @@ import test_csroot2caldb
 import test_cssens
 import test_csspec
 import test_cstsdist
+import test_cstsmapsplit
 import test_cstsmapmerge
 import test_csworkflow
 import test_csiactcopy
@@ -111,6 +112,7 @@ def test(installed=False, debug=False):
     suite_cssens       = test_cssens.Test()
     suite_csspec       = test_csspec.Test()
     suite_cstsdist     = test_cstsdist.Test()
+    suite_cstsmapsplit = test_cstsmapsplit.Test()
     suite_cstsmapmerge = test_cstsmapmerge.Test()
     suite_csworkflow   = test_csworkflow.Test()
  
@@ -128,6 +130,7 @@ def test(installed=False, debug=False):
     suite_cssens.set()
     suite_csspec.set()
     suite_cstsdist.set()
+    suite_cstsmapsplit.set()
     suite_cstsmapmerge.set()
     suite_csworkflow.set()
  
@@ -145,6 +148,7 @@ def test(installed=False, debug=False):
     suites.append(suite_cssens)
     suites.append(suite_csspec)
     suites.append(suite_cstsdist)
+    suites.append(suite_cstsmapsplit)
     suites.append(suite_cstsmapmerge)
     suites.append(suite_csworkflow)