Fabio Iocco / ctools

Commit 57c9b60cd61eb3fffe875ec4893f819aa2d25a23

Authored by Michael Mayer
1 parent 95b202bf

Added csmodelinfo and removed comment in csobsinfo

Inline Side-by-side
Showing 2 changed files with 511 additions and 1 deletions
cscripts/csmodelinfo.py 0 → 100644
  View file @57c9b60
  1 +#!/usr/bin/env python
  2 +# ==========================================================================
  3 +# Dump information about models into log file
  4 +#
  5 +# Copyright (C) 2015 Michael Mayer
  6 +#
  7 +# This program is free software: you can redistribute it and/or modify
  8 +# it under the terms of the GNU General Public License as published by
  9 +# the Free Software Foundation, either version 3 of the License, or
  10 +# (at your option) any later version.
  11 +#
  12 +# This program is distributed in the hope that it will be useful,
  13 +# but WITHOUT ANY WARRANTY; without even the implied warranty of
  14 +# MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
  15 +# GNU General Public License for more details.
  16 +#
  17 +# You should have received a copy of the GNU General Public License
  18 +# along with this program. If not, see <http://www.gnu.org/licenses/>.
  19 +#
  20 +# ==========================================================================
  21 +import gammalib
  22 +import ctools
  23 +import sys
  24 +
  25 +
  26 +# =============== #
  27 +# cshessobs class #
  28 +# =============== #
  29 +class csmodelinfo(ctools.cscript):
  30 + """
  31 + This class dumps information about a model container into a logfile or on screen.
  32 + This might be helpful for quick access to the model container, e.g. number of sources,
  33 + background models, free and fixed parameters etc.
  34 + The tool can also dump the models into a region file to display e.g. with ds9.
  35 + """
  36 + def __init__(self, *argv):
  37 + """
  38 + Constructor.
  39 + """
  40 + # Set name
  41 + self.name = "csmodelinfo"
  42 + self.version = "0.1.0"
  43 +
  44 + # Make sure that parfile exists
  45 + file = self.parfile()
  46 +
  47 + # Initialise application
  48 + if len(argv) == 0:
  49 + ctools.cscript.__init__(self, self.name, self.version)
  50 + elif len(argv) == 1:
  51 + ctools.cscript.__init__(self, self.name, self.version, *argv)
  52 + else:
  53 + raise TypeError("Invalid number of arguments given.")
  54 +
  55 + # Set logger properties
  56 + self.log_header()
  57 + self.log.date(True)
  58 +
  59 + # Return
  60 + return
  61 +
  62 + def __del__(self):
  63 + """
  64 + Destructor.
  65 + """
  66 + # Return
  67 + return
  68 +
  69 + def parfile(self):
  70 + """
  71 + Check if parfile exists. If parfile does not exist then create a
  72 + default parfile. This kluge avoids shipping the cscript with a parfile.
  73 + """
  74 +
  75 + # Set parfile name
  76 + parfile = self.name+".par"
  77 +
  78 + try:
  79 + pars = gammalib.GApplicationPars(parfile)
  80 + except:
  81 + # Signal that parfile was not found
  82 + print("Parfile \""+parfile+"\" not found. Create default parfile.")
  83 +
  84 + # Create default parfile
  85 + pars = gammalib.GApplicationPars()
  86 + pars.append(gammalib.GApplicationPar("inmodel","f","a","model.xml","","","Input model XML file"))
  87 + pars.append(gammalib.GApplicationPar("ds9file","f","h","NONE","","","DS9 region file containing models"))
  88 + pars.append(gammalib.GApplicationPar("pnt_type","s","h","cross","circle|box|diamond|cross|x|arrow|boxcircle","","Marker type for point sources"))
  89 + pars.append(gammalib.GApplicationPar("pnt_mark_size","i","h","12","","","Marker size for point sources"))
  90 + pars.append(gammalib.GApplicationPar("show_labels","b","h","yes","yes|no","","Add source labels?"))
  91 + pars.append(gammalib.GApplicationPar("width","i","h","2","","","Line width for regions"))
  92 + pars.append(gammalib.GApplicationPar("font","s","h","helvetica","","","Font for source labels"))
  93 + pars.append(gammalib.GApplicationPar("fontsize","i","h","12","","","Font size for source labels"))
  94 + pars.append(gammalib.GApplicationPar("fonttype","s","h","normal","normal|bold","","Normal or bold font?"))
  95 + pars.append(gammalib.GApplicationPar("fonttype2","s","h","roman","roman|italic","","Use roman or italic font?"))
  96 + pars.append(gammalib.GApplicationPar("show_ext_type","b","h","yes","yes|no","","Show type of extended model in source name"))
  97 + pars.append(gammalib.GApplicationPar("free_color","s","h","green","","","Color for sources with free parameters (any ds9 color or hex code)"))
  98 + pars.append(gammalib.GApplicationPar("fixed_color","s","h","magenta","","","Color for source with fixed parameters (any ds9 color or hex code)"))
  99 + pars.append_standard()
  100 + pars.append(gammalib.GApplicationPar("logfile","f","h","csmodelinfo.log","","","Log filename"))
  101 + pars.save(parfile)
  102 +
  103 + # Return
  104 + return
  105 +
  106 + def get_parameters(self):
  107 + """
  108 + Get parameters from parfile and setup the observation.
  109 + """
  110 + # Get models
  111 + self.models = gammalib.GModels(self["inmodel"].filename())
  112 +
  113 + # Get ds9 file
  114 + self.ds9file = self["ds9file"].filename()
  115 +
  116 + # Get hidden parameters for region file
  117 + self.pnt_type = self["pnt_type"].string()
  118 + self.pnt_mark_size = self["pnt_mark_size"].integer()
  119 + self.show_labels = self["show_labels"].boolean()
  120 + self.width = self["width"].integer()
  121 + self.font = self["font"].string()
  122 + self.fontsize = self["fontsize"].integer()
  123 + self.fonttype = self["fonttype"].string()
  124 + self.fonttype2 = self["fonttype2"].string()
  125 + self.show_ext_type = self["show_ext_type"].boolean()
  126 + self.free_color = self["free_color"].string()
  127 + self.fixed_color = self["fixed_color"].string()
  128 +
  129 + # Return
  130 + return
  131 +
  132 +
  133 + def run(self):
  134 + """
  135 + Run the script.
  136 + """
  137 + # Switch screen logging on in debug mode
  138 + if self.logDebug():
  139 + self.log.cout(True)
  140 +
  141 + # Get parameters
  142 + self.get_parameters()
  143 +
  144 + # Write input parameters into logger
  145 + if self.logTerse():
  146 + self.log_parameters()
  147 + self.log("\n")
  148 +
  149 + if self.logVerbose():
  150 + self.log("\n")
  151 + self.log.header1("Models in detail")
  152 + self.log("\n")
  153 + self.log(str(self.models))
  154 + self.log("\n")
  155 +
  156 + # Initialise output values
  157 + types = {}
  158 + instruments = {}
  159 + ts = {}
  160 + pars_at_limit = {}
  161 + free_src_pars = {}
  162 + n_par_total = 0
  163 + n_par_free = 0
  164 + n_par_fixed = 0
  165 + n_par_at_limit = 0
  166 + n_par_free_bkg = 0
  167 + n_par_free_src = 0
  168 + n_par_free_spec = 0
  169 + n_par_free_spat = 0
  170 + n_par_free_temp = 0
  171 +
  172 + # Loop over models and initialise dictionaries with model names
  173 + for model in self.models:
  174 + types[model.type()] = 0
  175 + instruments[model.instruments()] = 0
  176 + pars_at_limit[model.name()] = []
  177 +
  178 + # Loop over models again and fill values to be dumped
  179 + for model in self.models:
  180 +
  181 + # Model types
  182 + types[model.type()] += 1
  183 +
  184 + # Model instruments
  185 + instruments[model.instruments()] += 1
  186 +
  187 + # TS
  188 + if model.has_ts():
  189 + ts[model.name()] = model.ts()
  190 +
  191 + # Skymodel?
  192 + skymodel = (model.classname() == "GModelSky")
  193 + if skymodel:
  194 + free_src_pars[model.name()] = []
  195 +
  196 + # Loop over model parameters
  197 + for par in model:
  198 +
  199 + # Increment total parameters
  200 + n_par_total += 1
  201 +
  202 + # Check if parameter is at limit
  203 + if par.has_range():
  204 + if par.value() == par.min() or par.value() == par.max():
  205 + pars_at_limit[model.name()].append(par.name())
  206 + n_par_at_limit += 1
  207 +
  208 + # Check if model parameter is free/fixed and
  209 + # increment respective variable
  210 + if par.is_free():
  211 + n_par_free += 1
  212 + if skymodel:
  213 + free_src_pars[model.name()].append(par)
  214 + n_par_free_src += 1
  215 + else:
  216 + n_par_free_bkg += 1
  217 + else:
  218 + n_par_fixed += 1
  219 +
  220 + # If sky model, divide free parameters into spatial, spectral and temporal
  221 + if skymodel:
  222 + for par in model.spectral():
  223 + if par.is_free():
  224 + n_par_free_spec += 1
  225 + for par in model.spatial():
  226 + if par.is_free():
  227 + n_par_free_spat += 1
  228 + for par in model.temporal():
  229 + if par.is_free():
  230 + n_par_free_temp += 1
  231 +
  232 + # Log summary
  233 + self.log("\n")
  234 + self.log.header1("Summary")
  235 + self.log("\n")
  236 +
  237 + # Log instruments
  238 + self.log.header3("Instrument specific models")
  239 + for inst, n_inst in instruments.iteritems():
  240 + if inst == "":
  241 + inst = "All"
  242 + self.log.parformat(inst)
  243 + self.log(str(n_inst))
  244 + self.log("\n")
  245 +
  246 + # Log model types
  247 + self.log("\n")
  248 + self.log.header3("Model types")
  249 + for modeltype, n_types in types.iteritems():
  250 + self.log.parformat(modeltype)
  251 + self.log(str(n_types))
  252 + self.log("\n")
  253 +
  254 + # Log parameter information
  255 + self.log("\n")
  256 + self.log.header2("Parameter information")
  257 + self.log.parformat("All parameters")
  258 + self.log(str(n_par_total))
  259 + self.log("\n")
  260 + self.log.parformat("Fixed parameters")
  261 + self.log(str(n_par_fixed))
  262 + self.log("\n")
  263 + self.log.parformat("Free parameters (total)")
  264 + self.log(str(n_par_free))
  265 + self.log("\n")
  266 + self.log.parformat("Free background parameters")
  267 + self.log(str(n_par_free_bkg))
  268 + self.log("\n")
  269 + self.log.parformat("Free source parameters")
  270 + self.log(str(n_par_free_src))
  271 + self.log("\n")
  272 + if n_par_free_spec > 0:
  273 + self.log.parformat("Free spectral parameters")
  274 + self.log(str(n_par_free_spec))
  275 + self.log("\n")
  276 + if n_par_free_spat > 0:
  277 + self.log.parformat("Free spatial parameters")
  278 + self.log(str(n_par_free_spat))
  279 + self.log("\n")
  280 + if n_par_free_temp > 0:
  281 + self.log.parformat("Free temporal parameters")
  282 + self.log(str(n_par_free_temp))
  283 + self.log("\n")
  284 + self.log.parformat("Parameters at limit")
  285 + self.log(str(n_par_at_limit))
  286 + self.log("\n")
  287 + self.log("\n")
  288 +
  289 + # Log parameters at limit (if any)
  290 + if n_par_at_limit > 0:
  291 + self.log.header3("Parameters at limit")
  292 + for source, parameter in pars_at_limit.iteritems():
  293 + if len(parameter):
  294 + for par in parameter:
  295 + self.log.parformat(source)
  296 + self.log(par)
  297 + self.log("\n")
  298 +
  299 + # Log free parameters if requested
  300 + if self.logExplicit() and len(free_src_pars):
  301 + self.log("\n")
  302 + self.log.header2("Free source parameters")
  303 + for source, parameter in free_src_pars.iteritems():
  304 + if len(parameter):
  305 + self.log.header3(source)
  306 + for par in parameter:
  307 + self.log(str(par))
  308 + self.log("\n")
  309 + self.log("\n")
  310 +
  311 + # Log TS values if available
  312 + if len(ts):
  313 + self.log.header3("Test statistics")
  314 + for source,tsvalue in ts.iteritems():
  315 + self.log.parformat(source)
  316 + self.log("%.2f"%tsvalue)
  317 + self.log("\n")
  318 +
  319 + # Return
  320 + return
  321 +
  322 + def save(self):
  323 + """
  324 + Save pointings to ds9 region file if required
  325 + """
  326 +
  327 + # Check if ds9 file is valid
  328 + if not self.ds9file == "NONE":
  329 +
  330 + self.log("\n")
  331 + self.log.header1("Save region file")
  332 + self.log("\n")
  333 +
  334 + # Open file
  335 + f = open(self.ds9file,"w")
  336 +
  337 + # Write coordinate system
  338 + f.write("fk5\n")
  339 +
  340 + # Loop over models
  341 + for model in self.models:
  342 +
  343 + # Continue only if point source or extended source model
  344 + if model.type() == "PointSource" or model.type() == "ExtendedSource":
  345 +
  346 + line = self.model2ds9string(model)
  347 + if len(line):
  348 + f.write(line+"\n")
  349 +
  350 + # Logging for diffuse components
  351 + elif model.type() == "DiffuseSource":
  352 + self.log("Skipping diffuse model "+model.name()+"\n")
  353 +
  354 + # Logging for background components
  355 + else:
  356 + if self.logExplicit():
  357 + self.log("Skipping background model "+model.name()+"\n")
  358 +
  359 + # Close file
  360 + f.close()
  361 +
  362 + # Return
  363 + return
  364 +
  365 + def model2ds9string(self, model):
  366 + """
  367 + Converts model instance to ds9 region string
  368 + """
  369 +
  370 + # Initialise return string
  371 + ds9string = ""
  372 +
  373 + # Determine color
  374 + color = self.fixed_color
  375 + for par in model:
  376 + if par.is_free():
  377 + color = self.free_color
  378 + break
  379 +
  380 + # Initialise flag ot continue
  381 + is_valid = True
  382 +
  383 + # Retrieve model position
  384 + try:
  385 + modelpos = model.spatial().dir()
  386 + except AttributeError:
  387 + self.log("Skipping source "+model.name()+": Could not determine sky position\n")
  388 + is_valid = False
  389 +
  390 + # Continue only if position was found
  391 + if is_valid:
  392 +
  393 + # Retrieve model type and name
  394 + modeltype = model.type()
  395 + modelname = model.name()
  396 +
  397 + # Handle point source
  398 + if modeltype == "PointSource":
  399 +
  400 + ds9string += "point("
  401 + ds9string += str(modelpos.ra_deg())
  402 + ds9string += ","
  403 + ds9string += str(modelpos.dec_deg())
  404 + ds9string += ") # "
  405 + ds9string += "point="+self.pnt_type+" "
  406 + ds9string += str(self.pnt_mark_size)+" "
  407 +
  408 + # Hanlde extended sources
  409 + elif modeltype == "ExtendedSource":
  410 +
  411 + # Retrieve spatial model
  412 + spatial = model.spatial()
  413 + classname = spatial.classname()
  414 +
  415 + # Handle radial sources
  416 + if "Radial" in classname:
  417 +
  418 + # Retrieve short name of model class (e.g. "Gauss" or "Disk)
  419 + shorttype = classname.split("Radial")[-1]
  420 +
  421 + # Handle Disk and Shell model
  422 + if classname == "GModelSpatialRadialDisk" or classname == "GModelSpatialRadialShell":
  423 + size = spatial.radius()
  424 +
  425 + # Handle Gauss Model
  426 + elif classname == "GModelSpatialRadialGauss":
  427 + size = spatial.sigma()
  428 +
  429 + # Skip if source is unknown
  430 + else:
  431 + self.log("Skipping source "+model.name()+": Spatial model \""+classname+"\" unkown\n")
  432 + is_valid = False
  433 +
  434 + # Append radial ds9 string
  435 + if is_valid:
  436 + ds9string += "circle("
  437 + ds9string += str(modelpos.ra_deg())
  438 + ds9string += ","
  439 + ds9string += str(modelpos.dec_deg())
  440 + ds9string += ","
  441 + ds9string += str(size*3600.)+"\") #"
  442 +
  443 + # Handle elliptical sources
  444 + elif "Elliptical" in classname:
  445 +
  446 + # Retrieve short name and source size
  447 + shorttype = classname.split("Elliptical")[-1]
  448 + size1 = spatial.semimajor()
  449 + size2 = spatial.semiminor()
  450 + angle = spatial.posangle()
  451 +
  452 + # Append elliptical ds9 string
  453 + ds9string += "ellipse("
  454 + ds9string += str(modelpos.ra_deg())
  455 + ds9string += ","
  456 + ds9string += str(modelpos.dec_deg())
  457 + ds9string += ","
  458 + ds9string += str(size1*3600.)+"\""
  459 + ds9string += ","
  460 + ds9string += str(size2*3600.)+"\""
  461 + ds9string += ","
  462 + ds9string += str(angle+90.0)+") # "
  463 +
  464 + # Skip if source is neither radial nor elliptical
  465 + else:
  466 + self.log("Skipping source "+model.name()+": Spatial model neither \"Radial\" nor \"Elliptical\"\n")
  467 +
  468 + # Add short name to modelname
  469 + if self.show_ext_type:
  470 + modelname +=" ("+shorttype+")"
  471 +
  472 + # Add user defined attributes
  473 + ds9string += "color="+color+" "
  474 + ds9string += "width="+str(self.width)+" "
  475 + ds9string += "font=\""+self.font+" "+str(self.fontsize)+" "+self.fonttype+" "+self.fonttype2+"\" "
  476 + ds9string += "text={"+modelname+"}"
  477 +
  478 + # Return string
  479 + return ds9string
  480 +
  481 +
  482 + def execute(self):
  483 + """
  484 + Execute the script.
  485 + """
  486 + # Run the script
  487 + self.run()
  488 +
  489 + # Save ds9 file if required
  490 + self.save()
  491 +
  492 + # Return
  493 + return
  494 +
  495 +
  496 +# ======================== #
  497 +# Main routine entry point #
  498 +# ======================== #
  499 +if __name__ == '__main__':
  500 + """
  501 + Generates caldb entry.
  502 + """
  503 + # Create instance of application
  504 + app = csmodelinfo(sys.argv)
  505 +
  506 + # Open logfile
  507 + app.logFileOpen()
  508 +
  509 + # Execute application
  510 + app.execute()
... ...
cscripts/csobsinfo.py
  View file @57c9b60
... ... @@ -99,10 +99,10 @@ class csobsinfo(ctools.cscript):
99 99 Get parameters from parfile and setup the observation.
100 100 """
101 101 # Get parameters
102   - #self.inobs = self["inobs"].filename()
103 102 self.require_inobs("csobsinfo::get_parameters")
104 103 self.obs = self.get_observations(False)
105 104  
  105 + # Initialise object position
106 106 self.obj_dir = gammalib.GSkyDir()
107 107  
108 108 self.m_offset = self["offset"].boolean()
... ...