Add csmodelselect script (#1879)

Jürgen Knödlseder
1 parent 2cbd1b20
Showing 11 changed files with 655 additions and 25 deletions
ChangeLog
NEWS
cscripts/Makefile.am
cscripts/__init__.py.in
cscripts/csmodelselect.par
cscripts/csmodelselect.py
doc/source/users/reference_manual/csmodelselect.rst
doc/source/users/reference_manual/cspull.rst
doc/source/users/reference_manual/index.rst
test/test_csmodelselect.py
test/test_python_cscripts.py
-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)
 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
 --------
@@ -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 \
@@ -29,6 +29,7 @@ __all__ = [
     "cslightcrv",
     "csmodelinfo",
     "csmodelmerge",
+    "csmodelselect",
     "csobs2caldb",
     "csobsdef",
     "csobsinfo",
@@ -54,28 +55,29 @@ if cur_version &gt; 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
+#############################################################################
+#                                                                           #
+#                       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"
+#!/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 <http://www.gnu.org/licenses/>.
+#
+# ==========================================================================
+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()
+.. _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
@@ -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").
@@ -44,6 +44,7 @@ cscripts
    cslightcrv --- Computes lightcurve <cslightcrv>
    csmodelinfo --- Shows model container content <csmodelinfo>
    csmodelmerge --- Merges several model containers into one file <csmodelmerge>
+   csmodelselect --- Select models from model definition file <csmodelselect>
    csobsdef --- Generates observation definition file <csobsdef>
    csobsinfo --- Shows observation container content <csobsinfo>
    cspull --- Generates pull distribution <cspull>
+#! /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 <http://www.gnu.org/licenses/>.
+#
+# ==========================================================================
+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
@@ -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(),