From 63661f50eaf9b47ab13ac11dd48cec9e5f091584 Mon Sep 17 00:00:00 2001 From: Jürgen Knödlseder Date: Fri, 29 Jul 2016 18:20:46 +0200 Subject: [PATCH] Test all residual map algorithms in csresmap (#1828) --- cscripts/csresmap.py | 274 +++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++----------------------------------------------------------------------------------------------------------------------------------------------------------- test/test_csresmap.py | 7 ++++--- 2 files changed, 123 insertions(+), 158 deletions(-) diff --git a/cscripts/csresmap.py b/cscripts/csresmap.py index 65d07d7..c1cca42 100755 --- a/cscripts/csresmap.py +++ b/cscripts/csresmap.py @@ -28,12 +28,7 @@ import ctools # ============== # class csresmap(ctools.cscript): """ - Generates a residual map. - - This class implements the creation of a residual map. It derives from - the ctools.cscript class which provides support for parameter files, - command line arguments, and logging. In that way the Python script - behaves just as a regular ctool. + Generates a residual map """ # Constructor @@ -42,42 +37,18 @@ class csresmap(ctools.cscript): Constructor. """ # Set name - self._name = "csresmap" - self._version = "1.1.0" + self._name = 'csresmap' + self._version = '1.1.0' # Initialise class members - self._algorithm = "SUB" self._resmap = None - self._modcube = "NONE" self._use_maps = False self._skip_binning = False - self._publish = False - self._outfile = "" - self._use_maps = False - self._skip_binning = False - self._edisp = False - self._xref = 83.63 - self._yref = 22.01 - self._emin = 0.1 - self._emax = 100.0 - self._enumbins = 0 - self._ebinalg = "LOG" - self._coordsys = "CEL" - self._proj = "CAR" - self._nxpix = 200 - self._nypix = 200 - self._binsz = 0.02 - self._publish = False - self._chatter = 2 - self._clobber = True - self._debug = False - self._log_clients = False - - # Initialise observation container from constructor arguments. + + # Initialise observation container from constructor arguments self._obs, argv = self._set_input_obs(argv) - # Initialise application by calling the appropriate class - # constructor. + # Initialise application by calling the appropriate class constructor self._init_cscript(argv) # Return @@ -87,25 +58,29 @@ class csresmap(ctools.cscript): # Private methods def _get_parameters(self): """ - Get parameters from parfile and setup the observation. + Get parameters from parfile and setup the observation """ # Initialise some flags self._use_maps = False self._skip_binning = False - # First check if the inobs parameter is a counts cube - if self._obs.size() == 0 and self["inobs"].filename() != "NONE": - filename = gammalib.GFilename(self["inobs"].filename()) + # If there are no observations in the observation container then check + # if the "inobs" parameter is a counts cube. If the "inobs" parameter + # represents a FITS file and if this FITS file is a counts cube then + # set self._skip_binning=True + if self._obs.size() == 0 and self['inobs'].filename() != 'NONE': + filename = gammalib.GFilename(self['inobs'].filename()) if filename.is_fits(): cta = gammalib.GCTAObservation() cta.load(filename) - if cta.eventtype() == "CountsCube": + if cta.eventtype() == 'CountsCube': self._skip_binning = True - # If we have a counts cube, then ask whether we also have a model + # If we have a counts cube, then ask whether we also have a model cube. + # If a model cube name is given then set self._use_maps=True if self._skip_binning: - self._modcube = self["modcube"].filename() - if self._modcube != "NONE": + modcube = self['modcube'].filename() + if modcube != 'NONE': self._use_maps = True # If not two maps are given, proceed to set up observation @@ -113,56 +88,51 @@ class csresmap(ctools.cscript): # Set observation if not done before if self._obs.size() == 0: - self._require_inobs("csresmap.get_parameters()") + self._require_inobs('csresmap.get_parameters()') self._obs = self._get_observations() - # Check if we have exactly one binned CTA observation - if self._obs.size() == 1: - - if self._obs[0].classname() == "GCTAObservation": - if self._obs[0].eventtype() == "CountsCube": - # Skip ctbin step later on - self._skip_binning = True + # If we have exactly one binned CTA observation then signal that + # the binning can be skipped + if self._obs.size() == 1 and \ + self._obs[0].classname() == 'GCTAObservation' and \ + self._obs[0].eventtype() == 'CountsCube': + self._skip_binning = True # Set models if we have none if self._obs.models().size() == 0: - self._obs.models(self["inmodel"].filename()) + self._obs.models(self['inmodel'].filename()) - # Skip query for spatial parameters if a binning is provided in the observation + # If no binning is provided in the observation then query now + # the counts cube binning parameters if not self._skip_binning: - # Read other parameters - self._xref = self["xref"].real() - self._yref = self["yref"].real() - self._emin = self["emin"].real() - self._emax = self["emax"].real() - self._enumbins = self["enumbins"].integer() - self._ebinalg = self["ebinalg"].string() - self._coordsys = self["coordsys"].string() - self._proj = self["proj"].string() - self._nxpix = self["nxpix"].integer() - self._nypix = self["nypix"].integer() - self._binsz = self["binsz"].real() + self['xref'].real() + self['yref'].real() + self['emin'].real() + self['emax'].real() + self['enumbins'].integer() + self['ebinalg'].string() + self['coordsys'].string() + self['proj'].string() + self['nxpix'].integer() + self['nypix'].integer() + self['binsz'].real() - # Read energy dispersion flag - self._edisp = self["edisp"].boolean() - - # Read algorithm - self._algorithm = self["algorithm"].string() - - # Read standard parameters - self._publish = self["publish"].boolean() - self._chatter = self["chatter"].integer() - self._clobber = self["clobber"].boolean() - self._debug = self["debug"].boolean() - - # Read ahead output parameters + # Query parameters + self['edisp'].boolean() + self['algorithm'].string() + self['publish'].boolean() + self['chatter'].integer() + self['clobber'].boolean() + self['debug'].boolean() + + # Query ahead output parameters if (self._read_ahead()): - self["outmap"].filename() + self['outmap'].filename() # Write input parameters into logger if self._logTerse(): self._log_parameters() - self._log("\n") + self._log('\n') # Return return @@ -171,7 +141,7 @@ class csresmap(ctools.cscript): # Public methods def run(self): """ - Run the script. + Run the script """ # Switch screen logging on in debug mode if self._logDebug(): @@ -182,46 +152,48 @@ class csresmap(ctools.cscript): # Write observation into logger if self._logTerse(): - self._log("\n") - self._log.header1("Observation") + self._log('\n') + self._log.header1(gammalib.number('Observation',len(self._obs))) self._log(str(self._obs)) - self._log("\n") + self._log('\n') - # Use input file directly if given + # If a counts and model cube are specified then load them as sky map if self._use_maps: - countmap = gammalib.GSkyMap(self["inobs"].filename()) - modelmap = gammalib.GSkyMap(self._modcube) + countmap = gammalib.GSkyMap(self['inobs'].filename()) + modelmap = gammalib.GSkyMap(self['modcube'].filename()) + # ... otherwise build a counts cube and model cube else: - # ... + # Do not build counts cube if we have already one in the observation + # container if self._skip_binning: - cta_counts_cube = gammalib.GCTAEventCube(self._obs[0].events().clone()) + cta_counts_cube = gammalib.GCTAEventCube(self._obs[0].events()) - # ... + # ... otherwise generate one now from the event list else: # Write header if self._logTerse(): - self._log("\n") - self._log.header1("Generate binned map (ctbin)") + self._log('\n') + self._log.header1('Generate binned map (ctbin)') # Create countsmap bin = ctools.ctbin(self._obs) - bin["nxpix"].integer(self._nxpix) - bin["nypix"].integer(self._nypix) - bin["proj"].string(self._proj) - bin["coordsys"].string(self._coordsys) - bin["xref"].real(self._xref) - bin["yref"].real(self._yref) - bin["enumbins"].integer(self._enumbins) - bin["ebinalg"].string(self._ebinalg) - bin["emin"].real(self._emin) - bin["emax"].real(self._emax) - bin["binsz"].real(self._binsz) - bin["chatter"].integer(self._chatter) - bin["clobber"].boolean(self._clobber) - bin["debug"].boolean(self._debug) + bin['nxpix'] = self['nxpix'].integer() + bin['nypix'] = self['nypix'].integer() + bin['proj'] = self['proj'].string() + bin['coordsys'] = self['coordsys'].string() + bin['xref'] = self['xref'].real() + bin['yref'] = self['yref'].real() + bin['enumbins'] = self['enumbins'].integer() + bin['ebinalg'] = self['ebinalg'].string() + bin['emin'] = self['emin'].real() + bin['emax'] = self['emax'].real() + bin['binsz'] = self['binsz'].real() + bin['chatter'] = self['chatter'].integer() + bin['clobber'] = self['clobber'].boolean() + bin['debug'] = self['debug'].boolean() bin.run() # Retrieve counts cube @@ -232,63 +204,51 @@ class csresmap(ctools.cscript): # Write header if self._logTerse(): - self._log("\n") - self._log.header1("Generate model map (ctmodel)") + self._log('\n') + self._log.header1('Generate model map (ctmodel)') # Create model map model = ctools.ctmodel(self._obs) model.cube(cta_counts_cube) - model["chatter"].integer(self._chatter) - model["clobber"].boolean(self._clobber) - model["debug"].boolean(self._debug) - model["edisp"].boolean(self._edisp) + model['chatter'] = self['chatter'].integer() + model['clobber'] = self['clobber'].boolean() + model['debug'] = self['debug'].boolean() + model['edisp'] = self['edisp'].boolean() model.run() # Get model map into GSkyMap object modelmap = model.cube().counts().copy() - # Store counts map as residual map. Note that we need a - # special construct here to avoid memory leaks. This seems - # to be a SWIG feature as SWIG creates a new object when - # calling bin.cube() - #residualmap = bin.cube().counts() + # Store counts map as residual map. Note that we need a special + # construct here to avoid memory leaks. This seems to be a SWIG feature + # as SWIG creates a new object when calling bin.cube() self._resmap = countmap.copy() self._resmap.stack_maps() modelmap.stack_maps() - # Continue calculations depending on given algorithm - if self._algorithm == "SUB": + # Get residual map algorithm type + algorithm = self['algorithm'].string() - # Subtract maps + # Subtract maps + if algorithm == 'SUB': self._resmap -= modelmap - elif self._algorithm == "SUBDIV": - - # Subtract and divide by model map + # Subtract and divide by model map + elif algorithm == 'SUBDIV': self._resmap -= modelmap self._resmap /= modelmap - #for pixel in modelmap: - # if pixel != 0.0: - # pixel = 1.0/pixel - #self._resmap *= modelmap - elif self._algorithm == "SUBDIVSQRT": - - # subtract and divide by sqrt of model map + # Subtract and divide by sqrt of model map + elif algorithm == 'SUBDIVSQRT': self._resmap -= modelmap self._resmap /= modelmap.sqrt() - #for pixel in modelmap: - # if pixel != 0.0: - # pixel = 1.0/math.sqrt(pixel) - #self._resmap *= modelmap + # Raise exception if algorithm is unknown else: - - # Raise error if algorithm is unkown - raise TypeError("Algorithm \""+self._algorithm+"\" not known") + raise TypeError('Algorithm "'+algorithm+'" not known') # Optionally publish map - if self._publish: + if self['publish'].boolean(): self.publish() # Return @@ -296,7 +256,7 @@ class csresmap(ctools.cscript): def execute(self): """ - Execute the script. + Execute the script """ # Open logfile self.logFileOpen() @@ -315,55 +275,59 @@ class csresmap(ctools.cscript): def save(self): """ - Save residual map. + Save residual map """ # Write header if self._logTerse(): - self._log("\n") - self._log.header1("Save residual map") + self._log('\n') + self._log.header1('Save residual map') # Get outmap parameter - outmap = self["outmap"].filename() + outmap = self['outmap'].filename() # Continue only filename and residual map are valid if self._resmap != None: # Log file name if self._logTerse(): - self._log(gammalib.parformat("Residual map file")) + self._log(gammalib.parformat('Residual map file')) self._log(outmap.url()) - self._log("\n") + self._log('\n') # Save residual map - self._resmap.save(outmap, self._clobber) + self._resmap.save(outmap, self['clobber'].boolean()) # Return return - def publish(self, name=""): + def publish(self, name=''): """ Publish residual map. - Kwargs: - name: Name of residual map. + Parameters + ---------- + name : str, optional + Name of residual map """ # Write header if self._logTerse(): - self._log("\n") - self._log.header1("Publish residual map") + self._log('\n') + self._log.header1('Publish residual map') # Continue only if residual map is valid if self._resmap != None: # Set default name is user name is empty if not name: - user_name = self.name + user_name = self._name else: user_name = name # Log map name if self._logTerse(): - self._log("Publish residual map \""+user_name+"\".\n") + self._log(gammalib.parformat('Residual map name')) + self._log(user_name) + self._log('\n') # Publish map self._resmap.publish(user_name) @@ -373,10 +337,10 @@ class csresmap(ctools.cscript): def models(self, models): """ - Set model. + Set model """ # Copy models - self._obs.models(models.clone()) + self._obs.models(models) # Return return diff --git a/test/test_csresmap.py b/test/test_csresmap.py index 1944d38..392e372 100755 --- a/test/test_csresmap.py +++ b/test/test_csresmap.py @@ -76,7 +76,7 @@ class Test(test): ' emin=0.1 emax=100.0 enumbins=20'+ \ ' nxpix=50 nypix=50 binsz=0.02'+ \ ' coordsys="CEL" proj="CAR" xref=83.63 yref=22.01'+ \ - ' algorithm="SUBDIV"'+ \ + ' algorithm="SUBDIVSQRT"'+ \ ' logfile="csresmap_cmd1.log" chatter=1' # Check if execution of wrong command fails @@ -98,7 +98,7 @@ class Test(test): ' emin=0.1 emax=100.0 enumbins=20'+ \ ' nxpix=50 nypix=50 binsz=0.02'+ \ ' coordsys="CEL" proj="CAR" xref=83.63 yref=22.01'+ \ - ' algorithm="SUBDIV"'+ \ + ' algorithm="SUBDIVSQRT"'+ \ ' logfile="csresmap_cmd2.log" chatter=2' # Check if execution failed @@ -130,7 +130,7 @@ class Test(test): resmap['proj'] = 'CAR' resmap['xref'] = 83.63 resmap['yref'] = 22.01 - resmap['algorithm'] = 'SUBDIV' + resmap['algorithm'] = 'SUB' resmap['logfile'] = 'csresmap_py1.log' resmap['chatter'] = 2 @@ -157,6 +157,7 @@ class Test(test): resmap['algorithm'] = 'SUBDIV' resmap['logfile'] = 'csresmap_py2.log' resmap['chatter'] = 3 + resmap['publish'] = True # Run csresmap script resmap.execute() -- libgit2 0.21.4