diff --git a/ChangeLog b/ChangeLog index 3fa23ea..520a344 100644 --- a/ChangeLog +++ b/ChangeLog @@ -1,8 +1,9 @@ -2016-11-21 +2016-11-22 * Version 1.2.0 released ======================== + Add csmodelselect script (#1879) Add startindex parameter to ctobssim (#1878) Fix energy boundary handline in cssens (#1873) Add outcovmat parameter to ctlike (#1863) diff --git a/NEWS b/NEWS index b74f336..a592f0d 100644 --- a/NEWS +++ b/NEWS @@ -1,6 +1,6 @@ New Features and Important Changes in ctools 1.2.0 -21 November 2016 +22 November 2016 Introduction @@ -219,6 +219,9 @@ results on a VO hub (#1871). Fix energy boundary handling in cssens that has been corrupted during code refactoring (#1873). +Added csmodelselect script that allows selecting models from a model definition +XML file (#1879). + Examples -------- diff --git a/cscripts/Makefile.am b/cscripts/Makefile.am index fb22b0e..4be21d4 100644 --- a/cscripts/Makefile.am +++ b/cscripts/Makefile.am @@ -26,6 +26,7 @@ pythonscripts = $(builddir)/__init__.py \ $(srcdir)/cslightcrv.py \ $(srcdir)/csmodelinfo.py \ $(srcdir)/csmodelmerge.py \ + $(srcdir)/csmodelselect.py \ $(srcdir)/csobs2caldb.py \ $(srcdir)/csobsdef.py \ $(srcdir)/csobsinfo.py \ @@ -53,6 +54,7 @@ parfiles = $(srcdir)/cscaldb.par \ $(srcdir)/cslightcrv.par \ $(srcdir)/csmodelinfo.par \ $(srcdir)/csmodelmerge.par \ + $(srcdir)/csmodelselect.par \ $(srcdir)/csobs2caldb.par \ $(srcdir)/csobsdef.par \ $(srcdir)/csobsinfo.par \ @@ -77,6 +79,7 @@ linkscripts = cscaldb \ cslightcrv \ csmodelinfo \ csmodelmerge \ + csmodelselect \ csobs2caldb \ csobsdef \ csobsinfo \ @@ -102,6 +105,7 @@ testscripts = $(srcdir)/tests/__init__.py \ $(top_srcdir)/test/test_cscaldb.py \ $(top_srcdir)/test/test_csmodelinfo.py \ $(top_srcdir)/test/test_csmodelmerge.py \ + $(top_srcdir)/test/test_csmodelselect.py \ $(top_srcdir)/test/test_csobs2caldb.py \ $(top_srcdir)/test/test_csobsdef.py \ $(top_srcdir)/test/test_csobsinfo.py \ diff --git a/cscripts/__init__.py.in b/cscripts/__init__.py.in index 152bb3d..864149b 100644 --- a/cscripts/__init__.py.in +++ b/cscripts/__init__.py.in @@ -29,6 +29,7 @@ __all__ = [ "cslightcrv", "csmodelinfo", "csmodelmerge", + "csmodelselect", "csobs2caldb", "csobsdef", "csobsinfo", @@ -54,28 +55,29 @@ if cur_version > req_version: "csiactobs"]) # Import -from cscripts.cscaldb import cscaldb -from cscripts.csinfo import csinfo -from cscripts.cslightcrv import cslightcrv -from cscripts.csmodelinfo import csmodelinfo -from cscripts.csmodelmerge import csmodelmerge -from cscripts.csobs2caldb import csobs2caldb -from cscripts.csobsdef import csobsdef -from cscripts.csobsinfo import csobsinfo -from cscripts.cspull import cspull -from cscripts.csresmap import csresmap -from cscripts.csroot2caldb import csroot2caldb -from cscripts.cssens import cssens -from cscripts.csspec import csspec -from cscripts.cstsdist import cstsdist -from cscripts.cstsmapmerge import cstsmapmerge -from cscripts.cstsmapsplit import cstsmapsplit -from cscripts.csviscube import csviscube -from cscripts.csworkflow import csworkflow -from cscripts import obsutils -from cscripts import modutils -from cscripts import calutils -from cscripts import ioutils +from cscripts.cscaldb import cscaldb +from cscripts.csinfo import csinfo +from cscripts.cslightcrv import cslightcrv +from cscripts.csmodelinfo import csmodelinfo +from cscripts.csmodelmerge import csmodelmerge +from cscripts.csmodelselect import csmodelselect +from cscripts.csobs2caldb import csobs2caldb +from cscripts.csobsdef import csobsdef +from cscripts.csobsinfo import csobsinfo +from cscripts.cspull import cspull +from cscripts.csresmap import csresmap +from cscripts.csroot2caldb import csroot2caldb +from cscripts.cssens import cssens +from cscripts.csspec import csspec +from cscripts.cstsdist import cstsdist +from cscripts.cstsmapmerge import cstsmapmerge +from cscripts.cstsmapsplit import cstsmapsplit +from cscripts.csviscube import csviscube +from cscripts.csworkflow import csworkflow +from cscripts import obsutils +from cscripts import modutils +from cscripts import calutils +from cscripts import ioutils if cur_version > req_version: from cscripts.csfindobs import csfindobs from cscripts.csiactcopy import csiactcopy diff --git a/cscripts/csmodelselect.par b/cscripts/csmodelselect.par new file mode 100644 index 0000000..b6efac9 --- /dev/null +++ b/cscripts/csmodelselect.par @@ -0,0 +1,41 @@ +############################################################################# +# # +# CTA Science Analysis Tools # +# # +# ------------------------------------------------------------------------- # +# # +# File: csmodelselect.par # +# # +# Author: Juergen Knoedlseder # +# jknodlseder@irap.omp.eu # +# IRAP, Toulouse, France # +# # +# Purpose: Parameter file for the csmodelselect script # +# # +############################################################################# +# +# Input and output parameters +#============================ +inobs, f, a, obs.xml,,, "Input event list or observation definition XML file" +inmodel, f, a, inmodel.xml,,, "Input model definition XML file" +outmodel, f, a, outmodel.xml,,, "Output model definition XML file" + +# +# Script parameters +#================== +roilimit, r, h, 4.5,,, "Maximum RoI radius (deg)" +roimargin, r, h, 0.1,,, "Radial margin to be added to RoIs (deg)" +ethres, r, h, 0.1,,, "Energy threshold for source flux selection (TeV)" +fluxlimit, r, h, 1.0e-12,,, "Minimum source flux for selection (ph/cm2/s)" +tslimit, r, h, 10.0,,, "Minimum Test Statistic for selection" +fit_pos, b, h, yes,,, "Fit source positions?" +fit_shape, b, h, no,,, "Fit source shapes?" + +# +# Standard parameters +#==================== +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, csmodelselect.log,,, "Log filename" diff --git a/cscripts/csmodelselect.py b/cscripts/csmodelselect.py new file mode 100755 index 0000000..e873140 --- /dev/null +++ b/cscripts/csmodelselect.py @@ -0,0 +1,325 @@ +#!/usr/bin/env python +# ========================================================================== +# Selects model from a model definition XML file +# +# Copyright (C) 2016 Juergen Knoedlseder +# +# 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 gammalib +import ctools + + +# =================== # +# csmodelselect class # +# =================== # +class csmodelselect(ctools.cscript): + """ + Selects model from a model container + """ + + # Constructor + def __init__(self, *argv): + """ + Constructor + """ + # Set name + self._name = 'csmodelselect' + self._version = '1.2.0' + + # Initialise observation container from constructor arguments + self._obs, argv = self._set_input_obs(argv) + + # Initialise class members + self._models = gammalib.GModels() + + # 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 + """ + # Query input parameters + self['inobs'].filename() + self['inmodel'].filename() + + # Query hidden parameters + self['roilimit'].real() + self['roimargin'].real() + self['ethres'].real() + self['fluxlimit'].real() + self['tslimit'].real() + self['fit_pos'].boolean() + self['fit_shape'].boolean() + + # Query ahead output model filename + if self._read_ahead(): + self['outmodel'].filename() + + # If there are no observations in container then get them from the + # parameter file + if self._obs.size() == 0: + self._obs = self._get_observations(False) + + # Get models + self._models = gammalib.GModels(self['inmodel'].filename()) + + # Write input parameters into logger + self._log_parameters(gammalib.TERSE) + + # Return + return + + def _select_model(self, model, obs): + """ + Select model + + Parameters + ---------- + model : `~gammalib.GModel` + Model + obs : `~gammalib.GObservation` + Observation container + + Returns + ------- + select : bool + Model selection flag + """ + # Get selection parameters + roilimit = self['roilimit'].real() + roimargin = self['roimargin'].real() + ethres = self['ethres'].real() + fluxlimit = self['fluxlimit'].real() + tslimit = self['tslimit'].real() + + # Initialise selection flag to True + select = True + msg = 'Select by default' + + # If the model has a spatial component then check of overlap + if hasattr(model, 'spatial'): + + # Determine the model bounding box + model_bounds = model.spatial().region() + + # Initialise with an unselect model + select = False + msg = 'Exclude since outside any RoI' + + # Loop over all observations and check whether the model falls + # into one of them + for o in obs: + + # Get RoI centre and radius. The RoI radius is limited by + # roilimit. A margin given by roimargin is added to the RoI + # radius. + obs_centre = o.roi().centre().dir() + obs_radius = o.roi().radius() + if obs_radius > roilimit: + obs_radius = roilimit + obs_radius += roimargin + + # Set circular sky region + obs_bounds = gammalib.GSkyRegionCircle(obs_centre, obs_radius) + + # If model overlaps with RoI then signal overlap + if obs_bounds.overlaps(model_bounds): + select = True + msg = 'Select due to overlap with at least one RoI' + break + + # If model is selected and if model has a TS value then apply TS + # selection + if select and model.has_ts(): + if model.ts() < tslimit: + select = False + msg = 'Exclude since below TS limit (TS=%.3f)' % model.ts() + + + # If model is selected and if model is a sky model then apply flux + # limit selection + if select and model.classname() == 'GModelSky': + emin = gammalib.GEnergy(ethres, 'TeV') + emax = gammalib.GEnergy(1000.0, 'TeV') + flux = model.spectral().flux(emin, emax) + if flux < fluxlimit: + select = False + msg = 'Exclude since below flux limit (Flux=%.3e ph/cm2/s)' % flux + + # Log model selection + self._log_value(gammalib.NORMAL, model.name(), msg) + + # Return selection flag + return select + + def _set_model_parameters(self, model): + """ + Set model parameters + + Parameters + ---------- + model : `~gammalib.GModel` + Model + + Returns + ------- + model : `~gammalib.GModel` + Model with parameter set + """ + # Get selection parameters + fit_pos = self['fit_pos'].boolean() + fit_shape = self['fit_shape'].boolean() + + # If the model has a spatial component then set the spatial parameters + if hasattr(model, 'spatial'): + + # Loop over all parameters + for par in model: + + # Handle position parameters + if par.name() == 'RA' or par.name() == 'DEC': + if fit_pos: + par.free() + else: + par.fix() + + # Handle shape parameters + elif par.name() == 'Radius' or \ + par.name() == 'Sigma' or \ + par.name() == 'Width' or \ + par.name() == 'PA' or \ + par.name() == 'MajorRadius' or \ + par.name() == 'MinorRadius': + if fit_shape: + par.free() + else: + par.fix() + + # Return model + return model + + + # Public methods + 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() + + # Initialise empty model container for selected models + models = gammalib.GModels() + + # Write input observation container into logger + self._log_observations(gammalib.NORMAL, self._obs, 'Input observation'); + + # Write input models into logger + self._log_models(gammalib.VERBOSE, self._models, 'Input model') + + # Write header + self._log_header1(gammalib.NORMAL, 'Model selection') + + # Loop over model components + for model in self._models: + + # If model should be selected then set model parameters and append + # model to the output container + if self._select_model(model, self._obs): + + # Set model parameters + model = self._set_model_parameters(model) + + # Append model + models.append(model) + + # Copy selected model into selected models + self._models = models + + # Write selected models into logger + self._log_models(gammalib.VERBOSE, self._models, 'Selected model') + + # Return + return + + def save(self): + """ + Save model definition XML file + """ + # Write header + self._log_header1(gammalib.TERSE, 'Save models') + + # Get output filename in case it was not read ahead + outmodel = self['outmodel'].filename() + + # If file exists and clobber flag is false then raise an exception + if outmodel.exists() and not self._clobber: + msg = ('Cannot save "'+outmodel.url()+'": File already exists. ' + 'Use parameter clobber=yes to allow overwriting of files.') + raise RuntimeError(msg) + + # Otherwise log filename and save file + else: + # Log filename + self._log_value(gammalib.NORMAL, 'Model definition XML file', + outmodel.url()) + + # Save models + self._models.save(outmodel) + + # Return + return + + def execute(self): + """ + Execute the script + """ + # Open logfile + self.logFileOpen() + + # Read ahead output parameters + self._read_ahead(True) + + # Run the script + self.run() + + # Save model file if required + self.save() + + # Return + return + + +# ======================== # +# Main routine entry point # +# ======================== # +if __name__ == '__main__': + + # Create instance of application + app = csmodelselect(sys.argv) + + # Execute application + app.execute() diff --git a/doc/source/users/reference_manual/csmodelselect.rst b/doc/source/users/reference_manual/csmodelselect.rst new file mode 100644 index 0000000..722f637 --- /dev/null +++ b/doc/source/users/reference_manual/csmodelselect.rst @@ -0,0 +1,86 @@ +.. _csmodelselect: + +csmodelselect +============= + +Selects models from a model definition XML file. + + +Synopsis +-------- + +This script selects all models from a model definiton XML file that spatially +overlap with the Region of Interests (RoIs) of a list of observations and writes +these models into a new model definition XML file. In addition to the spatial +selection, the tool allows selecting only models with a photon flux and/or a +Test Statistic value above a given limit. The script also allows to specify +whether the source positions and/or source shapes should be fit. + + +General parameters +------------------ + +``inobs [file]`` + Input event list or observation definition XML file + +``inmodel [file]`` + Input model definition XML file + +``outmodel [file]`` + Output model definition XML file + +``(roilimit = 4.5) [real]`` + Maximum RoI radius (deg) + +``(roimargin = 0.1) [real]`` + Radial margin to be added to RoIs (deg) + +``(ethres = 0.1) [real]`` + Energy threshold for source flux selection (TeV) + +``(fluxlimit = 1.0e-12) [real]`` + Minimum source flux for selection (ph/cm2/s) + +``(tslimit = 10.0) [real]`` + Minimum Test Statistic for selection + +``(fit_pos = yes) [boolean]`` + Fit source positions? + +``(fit_shape = no) [boolean]`` + Fit source shapes? + + +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 model 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 = csmodelselect.log) [filename]`` + Log filename. + + +Related tools or scripts +------------------------ + +None diff --git a/doc/source/users/reference_manual/cspull.rst b/doc/source/users/reference_manual/cspull.rst index 8c45426..7be3549 100644 --- a/doc/source/users/reference_manual/cspull.rst +++ b/doc/source/users/reference_manual/cspull.rst @@ -124,7 +124,8 @@ Standard parameters Specifies whether an existing output file should be overwritten. ``(debug = no) [boolean]`` - Enables debug mode. In debug mode the executable will dump any log file output to the console. + 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"). diff --git a/doc/source/users/reference_manual/index.rst b/doc/source/users/reference_manual/index.rst index d83c1c4..073351e 100644 --- a/doc/source/users/reference_manual/index.rst +++ b/doc/source/users/reference_manual/index.rst @@ -44,6 +44,7 @@ cscripts cslightcrv --- Computes lightcurve csmodelinfo --- Shows model container content csmodelmerge --- Merges several model containers into one file + csmodelselect --- Select models from model definition file csobsdef --- Generates observation definition file csobsinfo --- Shows observation container content cspull --- Generates pull distribution diff --git a/test/test_csmodelselect.py b/test/test_csmodelselect.py new file mode 100755 index 0000000..2888214 --- /dev/null +++ b/test/test_csmodelselect.py @@ -0,0 +1,164 @@ +#! /usr/bin/env python +# ========================================================================== +# This scripts performs unit tests for the csmodelselect script +# +# Copyright (C) 2016 Juergen Knoedlseder +# +# 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 +from testing import test + + +# =================================== # +# Test class for csmodelselect script # +# =================================== # +class Test(test): + """ + Test class for csmodelselect script + + This test class makes unit tests for the csmodelselect script by using it + from the command line and from Python. + """ + + # Constructor + def __init__(self): + """ + Constructor. + """ + # Call base class constructor + test.__init__(self) + + # Return + return + + # Set test functions + def set(self): + """ + Set all test functions. + """ + # Set test name + self.name('csmodelselect') + + # Append tests + self.append(self._test_cmd, 'Test csmodelselect on command line') + self.append(self._test_python, 'Test csmodelselect from Python') + + # Return + return + + # Test csmodelselect on command line + def _test_cmd(self): + """ + Test csmodelselect on the command line + """ + # Set script name + csmodelselect = self._script('csmodelselect') + + # Setup csmodelselect command + cmd = csmodelselect+' inobs="'+self._events+'"'+ \ + ' inmodel="'+self._model+'"'+ \ + ' outmodel="csmodelselect_cmd1.xml"'+ \ + ' logfile="csmodelselect_cmd1.log" chatter=1' + + # Check if execution of wrong command fails + self.test_assert(self._execute('command_that_does_not_exist') != 0, + 'Self test of test script') + + # Check if execution was successful + self.test_assert(self._execute(cmd) == 0, + 'Check successful execution from command line') + + # Check model file + self._check_model_file('csmodelselect_cmd1.xml', 2) + + # Setup csmodelselect command + cmd = csmodelselect+' inobs="'+self._events+'"'+ \ + ' inmodel="model_that_does_not_exist.xml"'+ \ + ' outmodel="csmodelselect_cmd2.xml"'+ \ + ' logfile="csmodelselect_cmd2.log"' + + # Check if execution failed + self.test_assert(self._execute(cmd) != 0, + 'Check invalid input file when executed from command line') + + # Return + return + + # Test csmodelselect from Python + def _test_python(self): + """ + Test csmodelselect from Python + """ + # Allocate empty csmodelselect script + modelselect = cscripts.csmodelselect() + + # Check that saving does not nothing + modelselect['outmodel'] = 'csmodelselect_py0.xml' + modelselect['logfile'] = 'csmodelselect_py0.log' + modelselect.logFileOpen() + modelselect.save() + + # Check than an empty model definition file has been created + self._check_model_file('csmodelselect_py0.xml', 0) + + # Check that clearing does not lead to an exception or segfault + modelselect.clear() + + # Now set csmodelselect parameters + modelselect['inobs'] = self._events + modelselect['inmodel'] = self._model + modelselect['outmodel'] = 'csmodelselect_py1.xml' + modelselect['logfile'] = 'csmodelselect_py1.log' + modelselect['chatter'] = 2 + + # Run csmodelselect script and save models + modelselect.logFileOpen() # Make sure we get a log file + modelselect.run() + modelselect.save() + + # Check model file + self._check_model_file('csmodelselect_py1.xml', 2) + + # Execute csmodelselect script again, now with a higher chatter level + # than before + modelselect['outmodel'] = 'csmodelselect_py2.xml' + modelselect['logfile'] = 'csmodelselect_py2.log' + modelselect['chatter'] = 3 + modelselect.logFileOpen() # Needed to get a new log file + modelselect.execute() + + # Check model file + self._check_model_file('csmodelselect_py2.xml', 2) + + # Return + return + + # Check model file + def _check_model_file(self, filename, number): + """ + Check model file + """ + # Open models + models = gammalib.GModels(filename) + + # Check number of models + self.test_value(models.size(), number, 'Check for '+str(number)+ + ' models in XML file') + + # Return + return diff --git a/test/test_python_cscripts.py b/test/test_python_cscripts.py index e457a87..b0baaad 100755 --- a/test/test_python_cscripts.py +++ b/test/test_python_cscripts.py @@ -30,6 +30,7 @@ import test_csinfo import test_cslightcrv import test_csmodelinfo import test_csmodelmerge +import test_csmodelselect import test_csobs2caldb import test_csobsdef import test_csobsinfo @@ -114,6 +115,7 @@ def test(installed=False, debug=False): test_cslightcrv.Test(), test_csmodelinfo.Test(), test_csmodelmerge.Test(), + test_csmodelselect.Test(), test_csobs2caldb.Test(), test_csobsdef.Test(), test_csobsinfo.Test(),