From 16f51af681e5497d317a1af9873768710f7f38aa Mon Sep 17 00:00:00 2001 From: Michael Mayer Date: Thu, 30 Jun 2016 16:32:56 +0200 Subject: [PATCH] Added cstssplit (code, docs, test) --- cscripts/Makefile.am | 4 ++++ cscripts/__init__.py.in | 2 ++ cscripts/cstssplit.par | 55 +++++++++++++++++++++++++++++++++++++++++++++++++++++++ cscripts/cstssplit.py | 307 +++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++ doc/source/users/reference_manual/cstssplit.rst | 126 ++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++ doc/source/users/reference_manual/index.rst | 2 ++ test/test_cstssplit.py | 184 ++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++ test/test_python_cscripts.py | 4 ++++ 8 files changed, 684 insertions(+), 0 deletions(-) create mode 100644 cscripts/cstssplit.par create mode 100755 cscripts/cstssplit.py create mode 100644 doc/source/users/reference_manual/cstssplit.rst create mode 100755 test/test_cstssplit.py diff --git a/cscripts/Makefile.am b/cscripts/Makefile.am index 68f349b..7aeb394 100644 --- a/cscripts/Makefile.am +++ b/cscripts/Makefile.am @@ -22,6 +22,7 @@ pythonscripts = __init__.py \ csiactcopy.py \ csiactdata.py \ csiactobs.py \ + cstssplit.py \ obsutils.py # Parfiles to be distributed in $(prefix)/syspfiles @@ -43,6 +44,7 @@ parfiles = cscaldb.par \ csfindobs.par \ csiactcopy.par \ csiactdata.par \ + cstssplit.par \ csiactobs.par # Scripts for which symbolic links will be created in bindir @@ -65,6 +67,7 @@ linkscripts = cscaldb \ csfindobs \ csiactcopy \ csiactdata \ + cstssplit \ csiactobs # Test scripts to be distributed in $(pyexecdir)/cscripts/tests @@ -88,6 +91,7 @@ testscripts = tests/__init__.py \ ../test/test_csfindobs.py \ ../test/test_csiactcopy.py \ ../test/test_csiactdata.py \ + ../test/test_cstssplit.py \ ../test/test_csiactobs.py # Test data to be distributed in $(pyexecdir)/cscripts/tests/data diff --git a/cscripts/__init__.py.in b/cscripts/__init__.py.in index 7d299d2..59aff7b 100644 --- a/cscripts/__init__.py.in +++ b/cscripts/__init__.py.in @@ -40,6 +40,7 @@ __all__ = [ "cstsdist", "cstsmapmerge", "csworkflow", + "cstssplit", "obsutils" ] if cur_version > req_version: @@ -65,6 +66,7 @@ from cscripts.csspec import csspec from cscripts.cstsdist import cstsdist from cscripts.cstsmapmerge import cstsmapmerge from cscripts.csworkflow import csworkflow +from cscripts.cstssplit import cstssplit from cscripts import obsutils if cur_version > req_version: from cscripts.csfindobs import csfindobs diff --git a/cscripts/cstssplit.par b/cscripts/cstssplit.par new file mode 100644 index 0000000..2308969 --- /dev/null +++ b/cscripts/cstssplit.par @@ -0,0 +1,55 @@ +############################################################################# +# # +# CTA Science Analysis Tools # +# # +# ------------------------------------------------------------------------- # +# # +# File: cstssplit.par # +# Version: 1.1.0 # +# # +# Author: Michael Mayer # +# michael.mayer@physik.hu-berlin.de # +# HU Berlin # +# # +# Purpose: Parameter file for the cttsmap 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, "cstssplit.log",,, "Log filename" diff --git a/cscripts/cstssplit.py b/cscripts/cstssplit.py new file mode 100755 index 0000000..c45c3de --- /dev/null +++ b/cscripts/cstssplit.py @@ -0,0 +1,307 @@ +#!/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 . +# +# ========================================================================== +import sys +import os +import math +import gammalib +import ctools + + +# ============== # +# cstssplit class # +# ============== # +class cstssplit(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 + self._name = 'cstssplit' + self._version = '1.1.0' + 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: + 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 + + # Complete command by binning parameters + sliced_command = base_command + ' binmin='+str(binmin)+' binmax='+str(binmax) + sliced_command += ' outmap='+outmap + ' logfile='+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__': + """ + Generates splitted TS map commands + """ + + # Create instance of application + app = cstssplit(sys.argv) + + # Execute application + app.execute() + +# diff --git a/doc/source/users/reference_manual/cstssplit.rst b/doc/source/users/reference_manual/cstssplit.rst new file mode 100644 index 0000000..cd250bc --- /dev/null +++ b/doc/source/users/reference_manual/cstssplit.rst @@ -0,0 +1,126 @@ +.. _cstssplit: + +cstssplit +======= + +Generate commands to split TS map computation + + +Synopsis +-------- + +This tool allows to generate a command sequence in order to compute a +Test Statistic (TS) map using e.g. several machines. The command sequence +is stored in an ASCII file. Each line contains an individual command. +Each command will call the tool :ref:`cttsmap` with the appropriate parameters. +To split the computation, each TS map that is computed has a suffix to the filename +according to the number of the command. The resulting individual TS maps can then be +merged via :ref:`cstsmapmerge`. + + +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. + +``(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 [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 [string]`` + Projection method. + +``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 'false', 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. + + +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 = cttsmap.log) [string]`` + Name of log file. + + +Related tools or scripts +------------------------ + +:ref:`cttsmap` +:ref:`cstsmapmerge` + + diff --git a/doc/source/users/reference_manual/index.rst b/doc/source/users/reference_manual/index.rst index c4c2108..269b1c9 100644 --- a/doc/source/users/reference_manual/index.rst +++ b/doc/source/users/reference_manual/index.rst @@ -50,7 +50,9 @@ cscripts cssens --- Computes CTA sensitivity csspec --- Computes spectral points cstsdist --- Generates TS distribution + cstssplit --- Creates commands to slice the ts map computations cstsmapmerge --- Merges slices from ts map computations + Scripts to manage an IACT database diff --git a/test/test_cstssplit.py b/test/test_cstssplit.py new file mode 100755 index 0000000..66ef1f8 --- /dev/null +++ b/test/test_cstssplit.py @@ -0,0 +1,184 @@ +#! /usr/bin/env python +# ========================================================================== +# This scripts performs unit tests for the cstssplit 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 . +# +# ========================================================================== +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('cstssplit') + + # Append tests + self.append(self._test_cmd, "Test cstssplit on command line") + #self.append(self._test_python, "Test cstssplit from Python") + + # Return + return + + # Test cstssplit on command line + def _test_cmd(self): + """ + Test cstssplit on the command line. + """ + # Kluge to set the command (installed version has no README file) + if os.path.isfile("README.md"): + cstssplit = "../cscripts/cstssplit.py" + else: + cstssplit = "cstssplit" + + # Setup cstssplit command + cmd = cstssplit+' inobs="selected_events.fits"'+ \ + ' inmodel="data/crab.xml"'+ \ + ' srcname="Crab"'+ \ + ' caldb="irf"'+ \ + ' irf="cta_dummy_irf"'+ \ + ' 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 cstssplit command + cmd = cstssplit+' inobs="selected_events.fits"'+ \ + ' inmodel="data/crab.xml"'+ \ + ' srcname="Does_not_exist"'+ \ + ' caldb="irf"'+ \ + ' irf="cta_dummy_irf"'+ \ + ' 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 cstssplit from Python + def _test_python(self): + """ + Test cstssplit from Python. + """ + # Set-up cstssplit + tssplit = cscripts.cstssplit() + tssplit['inobs'] = 'obs_selected.xml' + tssplit['inmodel'] = 'data/crab.xml' + tssplit['xref'] = 83.6331 + tssplit['yref'] = 22.01 + tssplit['coordsys'] = 'CEL' + tssplit['proj'] = 'CAR' + tssplit['binsz'] = 0.05 + tssplit['nxpix'] = 5 + tssplit['nypix'] = 5 + tssplit['compute_null'] = False + tssplit['bins_per_job'] = 5 + tssplit['outfile'] = 'tsmap_cmd_py1.txt' + + + # Run cstssplit script and save run list + tssplit.logFileOpen() # Make sure we get a log file + tssplit.run() + tssplit.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(), + "Output of cstssplit exists") + + # Return + return diff --git a/test/test_python_cscripts.py b/test/test_python_cscripts.py index bd8d68c..231ccc4 100755 --- a/test/test_python_cscripts.py +++ b/test/test_python_cscripts.py @@ -42,6 +42,7 @@ import test_csworkflow import test_csiactcopy import test_csiactdata import test_csiactobs +import test_cstssplit # ================== # @@ -113,6 +114,7 @@ def test(installed=False, debug=False): suite_cstsdist = test_cstsdist.Test() suite_cstsmapmerge = test_cstsmapmerge.Test() suite_csworkflow = test_csworkflow.Test() + suite_cstssplit = test_cstssplit.Test() # Setup unit tests suite_cscaldb.set() @@ -130,6 +132,7 @@ def test(installed=False, debug=False): suite_cstsdist.set() suite_cstsmapmerge.set() suite_csworkflow.set() + suite_cstssplit.set() # Append tests to container suites.append(suite_cscaldb) @@ -147,6 +150,7 @@ def test(installed=False, debug=False): suites.append(suite_cstsdist) suites.append(suite_cstsmapmerge) suites.append(suite_csworkflow) + suites.append(suite_cstssplit) # Append tests for Python 2.6+ (the IACT cscripts depend on the json # module which is only available since Python 2.6+ -- libgit2 0.21.4