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

from scipy import optimize

import sys

import math

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

# csupperlimit class #

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

class csupperlimit(gl.GApplication):

    """

    This class implements an the computation of an upper limit value

    using a standard bayesian approach. 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    = "csupperlimit"

        self.version = "0.0.1"

        # Initialise some members

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

            self.obs = argv[0]

        else:      

            self.obs      = gl.GObservations()

        self.source = "CrabNebula"

        self.parname = "Prefactor"

        self.confidence = 0.95

        self.bestloglike = 0.0

        # 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("caldb","s","h","","","","Calibration database"))

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

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

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

            pars.append_standard()

            pars.save(parfile)

        # Return

        return

    def _func(self,value): 

        # get spectral model

        spec = self.obs.models()[self.source].spectral()

        # get parmeter value

        val = value[0]*spec[self.parname].scale()

        # ensure value to be above minimum boundary

        if val <= spec[self.parname].min():

            return 0.0-self.bestloglike-self.dloglike

        # set value

        spec[self.parname].value(val)

        # Recalculate likelihood

        like = ct.ctlike(self.obs)

        like.run()

        logl = like.opt().value()   

        # Return likelihood difference

        return (logl-self.bestloglike-self.dloglike)

    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 delta log-like

        # this has to be improved and made more flexible    

        if self.confidence == 0.95:

            self.dloglike = 3.84/2.0

    def execute(self):

        """

        Execute the script.

        """

        # Run the script

        self.run()

        # Return

        return

    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")

        # Check if observation container is already optimised

        # might add a flag in the future for that

        if self.obs.function().value() == 0.0:

            print "Optimising observation container"

            like = ct.ctlike(self.obs)

            like.run()

            self.bestloglike = like.opt().value()

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

        # Fix all parameters

        for model in self.obs.models():

            for par in model:

                par.fix()

        # get spectral parameters and set start value

        spectral = self.obs.models()[self.source].spectral()

        value = spectral[self.parname].value()

        error = spectral[self.parname].error()

        scale = spectral[self.parname].scale()

        # set start value to 2 sigma above fit value

        self.startpar = value / scale + 2* error/scale

        # calculate best loglike value if not done above

        if not self.bestloglike == 0.0: 

            like = ct.ctlike(self.obs)

            like.run()

            self.bestloglike = like.opt().value()

        # use scipy to find upper limit solution

        # The method used here is to determine the parameter value

        # where the likelihood function has changed by dloglike

        sol = optimize.root(self._func, self.startpar, method='hybr')

        # Test for convergence

        if sol.success:

            self.ulimit = sol.x[0]*spectral[self.parname].scale()

            if self.m_chatter>0:

                print "Upper limit of parameter \""+self.parname +"\" of source \""+self.source+"\" determined to "+str(self.ulimit)

        else:

            self.ulimit = 0.0 

            if self.m_chatter >0:

                print "Upper limit could not be determined"

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

# Main routine entry point #

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

if __name__ == '__main__':

    """

    Generates upperlimit for a source model.

    """

    # Create instance of application

    app = csupperlimit(sys.argv)

    # Open logfile

    app.logFileOpen()

    # Execute application

    app.execute()