ctools

#! /usr/bin/env python

# ==========================================================================

# This script generates the pull distribution for all model parameters.

#

# Copyright (C) 2011-2014 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 gammalib as gl

import ctools as ct

import sys

import math

# ============ #

# cslike class #

# ============ #

class cslike(gl.GApplication):

    """

    This class implements an auxiliary python interface to ctlike.

    It offers some additional functionality, like computing TS values,

    and storing covariance matrix. It derives from the 

    GammaLib::GApplication class which provides support for parameter

    files, command line arguments, and logging. In that way the Python

    script behaves just as a regular ctool.

    """

    def __init__(self, *argv):

        """

        Constructor.

        """

        # Set name

        self.name    = "cslike"

        self.version = "0.0.1"

        # Initialise some members

        if isinstance(argv[0],gl.GObservations):

            self.obs = argv[0]

        else:      

            self.obs      = gl.GObservations()

        self.outfile = ""

        self.m_value=0

        self.m_covar = None

        self.m_status = -1

        self.m_ts_method = "BIASED"

        self.ts = {}

        self.LogL0 = {}

        # Make sure that parfile exists

        file = self.parfile()

        # Initialise application

        if len(argv) == 1:

            gl.GApplication.__init__(self, self.name, self.version)

        elif len(argv) ==2:

            gl.GApplication.__init__(self, self.name, self.version, argv[1:])

        else:

            raise TypeError("Invalid number of arguments given.")

        # Set logger properties

        self.log_header()

        self.log.date(True)

        # Return

        return

    def __del__(self):

        """

        Destructor.

        """

        #  Write separator into logger

        if self.logTerse():

            self.log("\n")

        # Return

        return

    def parfile(self):

        """

        Check if parfile exists. If parfile does not exist then create a

        default parfile. This kluge avoids shipping the cscript with a parfile.

        """

        # Set parfile name

        parfile = self.name+".par"

        try:

            pars = gl.GApplicationPars(parfile)

        except:

            # Signal if parfile was not found

            sys.stdout.write("Parfile "+parfile+" not found. Create default parfile.\n")

            # Create default parfile

            pars = gl.GApplicationPars()

            pars.append(gl.GApplicationPar("infile", "f", "a", "events.fits","","","observation definition file"))

            pars.append(gl.GApplicationPar("stat","s","h","POISSON","","","Optimization statistics"))

            pars.append(gl.GApplicationPar("refit",  "b", "h", "no","","", "Do refitting?"))

            pars.append(gl.GApplicationPar("srcmdl","f","a","$CTOOLS/share/models/crab.xml","","","Source model"))

            pars.append(gl.GApplicationPar("outfile","f","a","crab_optmodel.xml","","","Output file name"))

            pars.append(gl.GApplicationPar("caldb","s","h","","","","Calibration database"))

            pars.append(gl.GApplicationPar("irf","s","a","cta_dummy_irf","","","Instrument response function"))

            pars.append(gl.GApplicationPar("TSmethod", "s", "h", "BIASED","BIASED|UNBIASED","","Method to compute TS values (BIASED is faster)"))

            pars.append(gl.GApplicationPar("edisp","b","h","no","","","Apply energy dispersion?"))

            pars.append(gl.GApplicationPar("logfile", "f", "h", "cslike.log","","", "Log filename"))

            pars.append_standard()

            pars.save(parfile)

        # Return

        return

    def _add_ts(self):

        """ appends TS value to xml file """ 

        xml = gl.GXml(self.m_outfile) 

        lib = xml.element("source_library", 0);        

        for src_name in self.ts.keys():    

            for i in range(lib.elements()):

                src_element = lib.element("source", i);

                if src_element.attribute("name") == src_name:

                    src_element.attribute("TS",str(self.ts[src_name]))

                    break

        xml.save(self.m_outfile)

        # append likelihood value to last line of the file

        f = open(self.m_outfile,"a")

        f.write("loglike="+str(self.m_value))

        f.close()

    def get_parameters(self):

        """

        Get parameters from parfile and setup the observation.

        """

        # Setup obs if it isn't already set

        if self.obs.size() == 0:

            infile = self["infile"].filename()

            try:      

                self.obs = gl.GObservations(infile)

            except:

                obs0 = gl.GCTAObservation(infile)

                self.obs.append(obs0)

            mdlfile = self["srcmdl"].filename()

            self.obs.models(mdlfile)

            self.m_caldb    = self["caldb"].string()

            self.m_irf      = self["irf"].string()

        else:

            self.m_caldb = ""

            self.m_irf = ""

        # get hidden parameters

        self.m_stat = self["stat"].string()

        self.m_refit = self["refit"].boolean()

        self.m_edisp = self["edisp"].boolean()

        self.m_chatter = self["chatter"].integer()

        self.m_debug = self["debug"].boolean()

        self.m_clobber = self["clobber"].boolean()

        self.m_mode = self["mode"].string()

        self.m_logfile = self["logfile"].filename()

        # get TS calculation method

        self.m_ts_method = self["TSmethod"].string()

        # get outfile

        if self.outfile == "":

            self.m_outfile  = self["outfile"].filename() 

        else:

            self.m_outfile = self.outfile 

    def execute(self):

        """

        Execute the script.

        """

        # Run the script

        self.run()

        # Return

        return

    def covar(self):

        """

        Return covariance matrix

        """

        return self.m_covar.clone()

    def obs(self):

        """

        Return observations

        """

        return self.obs.clone()    

    def run(self):

        """

        Run the fit.

        """

        # Switch screen logging on in debug mode

        if self.logDebug():

            self.log.cout(True)

        # Get parameters

        self.get_parameters()

        #  Write input parameters into logger

        if self.logTerse():

            self.log_parameters()

            self.log("\n")

        # Write observation into logger

        if self.logTerse():

            self.log("\n")

            self.log.header1("Observation")

            self.log(str(self.obs))

            self.log("\n")

        # Write header

        if self.logTerse():

            self.log("\n")

            self.log.header1("Generate pull distribution")

        # create ctlike instance with given parameters

        like = ct.ctlike(self.obs)

        like["stat"].string(self.m_stat)

        like["caldb"].string(self.m_caldb)

        like["irf"].string(self.m_irf)

        like["clobber"].boolean(self.m_clobber)

        like["debug"].boolean(self.m_debug)

        like["mode"].string(self.m_mode)

        like["logfile"].filename(self.m_logfile)

        like["edisp"].boolean(self.m_edisp)

        like["chatter"].integer(self["chatter"].integer())

        # Perform the fit

        if self["chatter"].integer()>0:

            print "Running ctlike ..."

        like.run()

        LogL1 = like.opt().value()

        # Store likelihood value and fit status

        self.m_value = LogL1

        self.m_status = like.opt().status()

        # Get covariance matrix if fit has converged

        try:

            self.m_covar = like.obs().function().curvature().invert().clone()

        except:

            print "WARNING: Covariance matrix not determined successfully, FitStatus = "+str(self.status),sys._getframe().f_code.co_name

        # save best fit results if outfile is given

        if not self.m_outfile == "":

            like.obs().models().save(self.m_outfile)

        # assign optimised models   

        self.obs.models(like.obs().models().clone())

        # Find sources with free parameters to calculate their TS values

        freemodelnames = []

        for m in like.obs().models():

            for par in m:

                if par.is_free():

                    freemodelnames.append(m.name())

                    break  

            # Fix parameters to speed up computations: Method is inspired

            # from Fermi Science tools, but might lead to biased results

            if self.m_ts_method == "BIASED":        

                for par in m:

                    par.fix()

        # save initial models for TS computation            

        bck_models = like.obs().models().clone()

        # Loop over models with free parameters

        for name in freemodelnames: 

            # exclude background from TS computation

            if isinstance(like.obs().models()[name],gl.GModelData):

                continue 

            # remove respective model from container                

            like.obs().models().remove(name)

            # rerun ctlike without model

            like0 = ct.ctlike(like.obs())

            like0.run()

            # get loglike value

            loglike0 = like0.opt().value()

            TS = 2* (loglike0-LogL1)

            # append TS value to dictionary

            self.ts[name]=TS

            # store likelihood value

            self.LogL0[name]=loglike0

            if self.m_chatter > 0:

                print "TS value of "+name+": "+str(TS)

            # Reset models to best-fit state

            like.obs().models(bck_models)

        # Add TS values to optimised model files   

        if not self.m_outfile == "":

            self._add_ts()

# ======================== #

# Main routine entry point #

# ======================== #

if __name__ == '__main__':

    """

    Runs ctlike with some extensions.

    """

    # Create instance of application

    app = cslike(sys.argv)

    # Open logfile

    app.logFileOpen()

    # Execute application

    app.execute()