Andreas Specovius / ctools

Commit 63661f50eaf9b47ab13ac11dd48cec9e5f091584

Authored by Jürgen Knödlseder
1 parent 49bfb301

Test all residual map algorithms in csresmap (#1828)

Inline Side-by-side
Showing 2 changed files with 123 additions and 158 deletions
cscripts/csresmap.py
  View file @63661f5
@@ -28,12 +28,7 @@ import ctools @@ -28,12 +28,7 @@ import ctools
28 # ============== # 28 # ============== #
29 class csresmap(ctools.cscript): 29 class csresmap(ctools.cscript):
30 """ 30 """
31 - Generates a residual map.  
32 -  
33 - This class implements the creation of a residual map. It derives from  
34 - the ctools.cscript class which provides support for parameter files,  
35 - command line arguments, and logging. In that way the Python script  
36 - behaves just as a regular ctool. 31 + Generates a residual map
37 """ 32 """
38 33
39 # Constructor 34 # Constructor
@@ -42,42 +37,18 @@ class csresmap(ctools.cscript): @@ -42,42 +37,18 @@ class csresmap(ctools.cscript):
42 Constructor. 37 Constructor.
43 """ 38 """
44 # Set name 39 # Set name
45 - self._name = "csresmap"  
46 - self._version = "1.1.0" 40 + self._name = 'csresmap'
  41 + self._version = '1.1.0'
47 42
48 # Initialise class members 43 # Initialise class members
49 - self._algorithm = "SUB"  
50 self._resmap = None 44 self._resmap = None
51 - self._modcube = "NONE"  
52 self._use_maps = False 45 self._use_maps = False
53 self._skip_binning = False 46 self._skip_binning = False
54 - self._publish = False  
55 - self._outfile = ""  
56 - self._use_maps = False  
57 - self._skip_binning = False  
58 - self._edisp = False  
59 - self._xref = 83.63  
60 - self._yref = 22.01  
61 - self._emin = 0.1  
62 - self._emax = 100.0  
63 - self._enumbins = 0  
64 - self._ebinalg = "LOG"  
65 - self._coordsys = "CEL"  
66 - self._proj = "CAR"  
67 - self._nxpix = 200  
68 - self._nypix = 200  
69 - self._binsz = 0.02  
70 - self._publish = False  
71 - self._chatter = 2  
72 - self._clobber = True  
73 - self._debug = False  
74 - self._log_clients = False  
75 -  
76 - # Initialise observation container from constructor arguments. 47 +
  48 + # Initialise observation container from constructor arguments
77 self._obs, argv = self._set_input_obs(argv) 49 self._obs, argv = self._set_input_obs(argv)
78 50
79 - # Initialise application by calling the appropriate class  
80 - # constructor. 51 + # Initialise application by calling the appropriate class constructor
81 self._init_cscript(argv) 52 self._init_cscript(argv)
82 53
83 # Return 54 # Return
@@ -87,25 +58,29 @@ class csresmap(ctools.cscript): @@ -87,25 +58,29 @@ class csresmap(ctools.cscript):
87 # Private methods 58 # Private methods
88 def _get_parameters(self): 59 def _get_parameters(self):
89 """ 60 """
90 - Get parameters from parfile and setup the observation. 61 + Get parameters from parfile and setup the observation
91 """ 62 """
92 # Initialise some flags 63 # Initialise some flags
93 self._use_maps = False 64 self._use_maps = False
94 self._skip_binning = False 65 self._skip_binning = False
95 66
96 - # First check if the inobs parameter is a counts cube  
97 - if self._obs.size() == 0 and self["inobs"].filename() != "NONE":  
98 - filename = gammalib.GFilename(self["inobs"].filename()) 67 + # If there are no observations in the observation container then check
  68 + # if the "inobs" parameter is a counts cube. If the "inobs" parameter
  69 + # represents a FITS file and if this FITS file is a counts cube then
  70 + # set self._skip_binning=True
  71 + if self._obs.size() == 0 and self['inobs'].filename() != 'NONE':
  72 + filename = gammalib.GFilename(self['inobs'].filename())
99 if filename.is_fits(): 73 if filename.is_fits():
100 cta = gammalib.GCTAObservation() 74 cta = gammalib.GCTAObservation()
101 cta.load(filename) 75 cta.load(filename)
102 - if cta.eventtype() == "CountsCube": 76 + if cta.eventtype() == 'CountsCube':
103 self._skip_binning = True 77 self._skip_binning = True
104 78
105 - # If we have a counts cube, then ask whether we also have a model 79 + # If we have a counts cube, then ask whether we also have a model cube.
  80 + # If a model cube name is given then set self._use_maps=True
106 if self._skip_binning: 81 if self._skip_binning:
107 - self._modcube = self["modcube"].filename()  
108 - if self._modcube != "NONE": 82 + modcube = self['modcube'].filename()
  83 + if modcube != 'NONE':
109 self._use_maps = True 84 self._use_maps = True
110 85
111 # If not two maps are given, proceed to set up observation 86 # If not two maps are given, proceed to set up observation
@@ -113,56 +88,51 @@ class csresmap(ctools.cscript): @@ -113,56 +88,51 @@ class csresmap(ctools.cscript):
113 88
114 # Set observation if not done before 89 # Set observation if not done before
115 if self._obs.size() == 0: 90 if self._obs.size() == 0:
116 - self._require_inobs("csresmap.get_parameters()") 91 + self._require_inobs('csresmap.get_parameters()')
117 self._obs = self._get_observations() 92 self._obs = self._get_observations()
118 93
119 - # Check if we have exactly one binned CTA observation  
120 - if self._obs.size() == 1:  
121 -  
122 - if self._obs[0].classname() == "GCTAObservation":  
123 - if self._obs[0].eventtype() == "CountsCube":  
124 - # Skip ctbin step later on  
125 - self._skip_binning = True 94 + # If we have exactly one binned CTA observation then signal that
  95 + # the binning can be skipped
  96 + if self._obs.size() == 1 and \
  97 + self._obs[0].classname() == 'GCTAObservation' and \
  98 + self._obs[0].eventtype() == 'CountsCube':
  99 + self._skip_binning = True
126 100
127 # Set models if we have none 101 # Set models if we have none
128 if self._obs.models().size() == 0: 102 if self._obs.models().size() == 0:
129 - self._obs.models(self["inmodel"].filename()) 103 + self._obs.models(self['inmodel'].filename())
130 104
131 - # Skip query for spatial parameters if a binning is provided in the observation 105 + # If no binning is provided in the observation then query now
  106 + # the counts cube binning parameters
132 if not self._skip_binning: 107 if not self._skip_binning:
133 - # Read other parameters  
134 - self._xref = self["xref"].real()  
135 - self._yref = self["yref"].real()  
136 - self._emin = self["emin"].real()  
137 - self._emax = self["emax"].real()  
138 - self._enumbins = self["enumbins"].integer()  
139 - self._ebinalg = self["ebinalg"].string()  
140 - self._coordsys = self["coordsys"].string()  
141 - self._proj = self["proj"].string()  
142 - self._nxpix = self["nxpix"].integer()  
143 - self._nypix = self["nypix"].integer()  
144 - self._binsz = self["binsz"].real() 108 + self['xref'].real()
  109 + self['yref'].real()
  110 + self['emin'].real()
  111 + self['emax'].real()
  112 + self['enumbins'].integer()
  113 + self['ebinalg'].string()
  114 + self['coordsys'].string()
  115 + self['proj'].string()
  116 + self['nxpix'].integer()
  117 + self['nypix'].integer()
  118 + self['binsz'].real()
145 119
146 - # Read energy dispersion flag  
147 - self._edisp = self["edisp"].boolean()  
148 -  
149 - # Read algorithm  
150 - self._algorithm = self["algorithm"].string()  
151 -  
152 - # Read standard parameters  
153 - self._publish = self["publish"].boolean()  
154 - self._chatter = self["chatter"].integer()  
155 - self._clobber = self["clobber"].boolean()  
156 - self._debug = self["debug"].boolean()  
157 -  
158 - # Read ahead output parameters 120 + # Query parameters
  121 + self['edisp'].boolean()
  122 + self['algorithm'].string()
  123 + self['publish'].boolean()
  124 + self['chatter'].integer()
  125 + self['clobber'].boolean()
  126 + self['debug'].boolean()
  127 +
  128 + # Query ahead output parameters
159 if (self._read_ahead()): 129 if (self._read_ahead()):
160 - self["outmap"].filename() 130 + self['outmap'].filename()
161 131
162 # Write input parameters into logger 132 # Write input parameters into logger
163 if self._logTerse(): 133 if self._logTerse():
164 self._log_parameters() 134 self._log_parameters()
165 - self._log("\n") 135 + self._log('\n')
166 136
167 # Return 137 # Return
168 return 138 return
@@ -171,7 +141,7 @@ class csresmap(ctools.cscript): @@ -171,7 +141,7 @@ class csresmap(ctools.cscript):
171 # Public methods 141 # Public methods
172 def run(self): 142 def run(self):
173 """ 143 """
174 - Run the script. 144 + Run the script
175 """ 145 """
176 # Switch screen logging on in debug mode 146 # Switch screen logging on in debug mode
177 if self._logDebug(): 147 if self._logDebug():
@@ -182,46 +152,48 @@ class csresmap(ctools.cscript): @@ -182,46 +152,48 @@ class csresmap(ctools.cscript):
182 152
183 # Write observation into logger 153 # Write observation into logger
184 if self._logTerse(): 154 if self._logTerse():
185 - self._log("\n")  
186 - self._log.header1("Observation") 155 + self._log('\n')
  156 + self._log.header1(gammalib.number('Observation',len(self._obs)))
187 self._log(str(self._obs)) 157 self._log(str(self._obs))
188 - self._log("\n") 158 + self._log('\n')
189 159
190 - # Use input file directly if given 160 + # If a counts and model cube are specified then load them as sky map
191 if self._use_maps: 161 if self._use_maps:
192 - countmap = gammalib.GSkyMap(self["inobs"].filename())  
193 - modelmap = gammalib.GSkyMap(self._modcube) 162 + countmap = gammalib.GSkyMap(self['inobs'].filename())
  163 + modelmap = gammalib.GSkyMap(self['modcube'].filename())
194 164
  165 + # ... otherwise build a counts cube and model cube
195 else: 166 else:
196 167
197 - # ... 168 + # Do not build counts cube if we have already one in the observation
  169 + # container
198 if self._skip_binning: 170 if self._skip_binning:
199 - cta_counts_cube = gammalib.GCTAEventCube(self._obs[0].events().clone()) 171 + cta_counts_cube = gammalib.GCTAEventCube(self._obs[0].events())
200 172
201 - # ... 173 + # ... otherwise generate one now from the event list
202 else: 174 else:
203 175
204 # Write header 176 # Write header
205 if self._logTerse(): 177 if self._logTerse():
206 - self._log("\n")  
207 - self._log.header1("Generate binned map (ctbin)") 178 + self._log('\n')
  179 + self._log.header1('Generate binned map (ctbin)')
208 180
209 # Create countsmap 181 # Create countsmap
210 bin = ctools.ctbin(self._obs) 182 bin = ctools.ctbin(self._obs)
211 - bin["nxpix"].integer(self._nxpix)  
212 - bin["nypix"].integer(self._nypix)  
213 - bin["proj"].string(self._proj)  
214 - bin["coordsys"].string(self._coordsys)  
215 - bin["xref"].real(self._xref)  
216 - bin["yref"].real(self._yref)  
217 - bin["enumbins"].integer(self._enumbins)  
218 - bin["ebinalg"].string(self._ebinalg)  
219 - bin["emin"].real(self._emin)  
220 - bin["emax"].real(self._emax)  
221 - bin["binsz"].real(self._binsz)  
222 - bin["chatter"].integer(self._chatter)  
223 - bin["clobber"].boolean(self._clobber)  
224 - bin["debug"].boolean(self._debug) 183 + bin['nxpix'] = self['nxpix'].integer()
  184 + bin['nypix'] = self['nypix'].integer()
  185 + bin['proj'] = self['proj'].string()
  186 + bin['coordsys'] = self['coordsys'].string()
  187 + bin['xref'] = self['xref'].real()
  188 + bin['yref'] = self['yref'].real()
  189 + bin['enumbins'] = self['enumbins'].integer()
  190 + bin['ebinalg'] = self['ebinalg'].string()
  191 + bin['emin'] = self['emin'].real()
  192 + bin['emax'] = self['emax'].real()
  193 + bin['binsz'] = self['binsz'].real()
  194 + bin['chatter'] = self['chatter'].integer()
  195 + bin['clobber'] = self['clobber'].boolean()
  196 + bin['debug'] = self['debug'].boolean()
225 bin.run() 197 bin.run()
226 198
227 # Retrieve counts cube 199 # Retrieve counts cube
@@ -232,63 +204,51 @@ class csresmap(ctools.cscript): @@ -232,63 +204,51 @@ class csresmap(ctools.cscript):
232 204
233 # Write header 205 # Write header
234 if self._logTerse(): 206 if self._logTerse():
235 - self._log("\n")  
236 - self._log.header1("Generate model map (ctmodel)") 207 + self._log('\n')
  208 + self._log.header1('Generate model map (ctmodel)')
237 209
238 # Create model map 210 # Create model map
239 model = ctools.ctmodel(self._obs) 211 model = ctools.ctmodel(self._obs)
240 model.cube(cta_counts_cube) 212 model.cube(cta_counts_cube)
241 - model["chatter"].integer(self._chatter)  
242 - model["clobber"].boolean(self._clobber)  
243 - model["debug"].boolean(self._debug)  
244 - model["edisp"].boolean(self._edisp) 213 + model['chatter'] = self['chatter'].integer()
  214 + model['clobber'] = self['clobber'].boolean()
  215 + model['debug'] = self['debug'].boolean()
  216 + model['edisp'] = self['edisp'].boolean()
245 model.run() 217 model.run()
246 218
247 # Get model map into GSkyMap object 219 # Get model map into GSkyMap object
248 modelmap = model.cube().counts().copy() 220 modelmap = model.cube().counts().copy()
249 221
250 - # Store counts map as residual map. Note that we need a  
251 - # special construct here to avoid memory leaks. This seems  
252 - # to be a SWIG feature as SWIG creates a new object when  
253 - # calling bin.cube()  
254 - #residualmap = bin.cube().counts() 222 + # Store counts map as residual map. Note that we need a special
  223 + # construct here to avoid memory leaks. This seems to be a SWIG feature
  224 + # as SWIG creates a new object when calling bin.cube()
255 self._resmap = countmap.copy() 225 self._resmap = countmap.copy()
256 self._resmap.stack_maps() 226 self._resmap.stack_maps()
257 modelmap.stack_maps() 227 modelmap.stack_maps()
258 228
259 - # Continue calculations depending on given algorithm  
260 - if self._algorithm == "SUB": 229 + # Get residual map algorithm type
  230 + algorithm = self['algorithm'].string()
261 231
262 - # Subtract maps 232 + # Subtract maps
  233 + if algorithm == 'SUB':
263 self._resmap -= modelmap 234 self._resmap -= modelmap
264 235
265 - elif self._algorithm == "SUBDIV":  
266 -  
267 - # Subtract and divide by model map 236 + # Subtract and divide by model map
  237 + elif algorithm == 'SUBDIV':
268 self._resmap -= modelmap 238 self._resmap -= modelmap
269 self._resmap /= modelmap 239 self._resmap /= modelmap
270 - #for pixel in modelmap:  
271 - # if pixel != 0.0:  
272 - # pixel = 1.0/pixel  
273 - #self._resmap *= modelmap  
274 240
275 - elif self._algorithm == "SUBDIVSQRT":  
276 -  
277 - # subtract and divide by sqrt of model map 241 + # Subtract and divide by sqrt of model map
  242 + elif algorithm == 'SUBDIVSQRT':
278 self._resmap -= modelmap 243 self._resmap -= modelmap
279 self._resmap /= modelmap.sqrt() 244 self._resmap /= modelmap.sqrt()
280 - #for pixel in modelmap:  
281 - # if pixel != 0.0:  
282 - # pixel = 1.0/math.sqrt(pixel)  
283 - #self._resmap *= modelmap  
284 245
  246 + # Raise exception if algorithm is unknown
285 else: 247 else:
286 -  
287 - # Raise error if algorithm is unkown  
288 - raise TypeError("Algorithm \""+self._algorithm+"\" not known") 248 + raise TypeError('Algorithm "'+algorithm+'" not known')
289 249
290 # Optionally publish map 250 # Optionally publish map
291 - if self._publish: 251 + if self['publish'].boolean():
292 self.publish() 252 self.publish()
293 253
294 # Return 254 # Return
@@ -296,7 +256,7 @@ class csresmap(ctools.cscript): @@ -296,7 +256,7 @@ class csresmap(ctools.cscript):
296 256
297 def execute(self): 257 def execute(self):
298 """ 258 """
299 - Execute the script. 259 + Execute the script
300 """ 260 """
301 # Open logfile 261 # Open logfile
302 self.logFileOpen() 262 self.logFileOpen()
@@ -315,55 +275,59 @@ class csresmap(ctools.cscript): @@ -315,55 +275,59 @@ class csresmap(ctools.cscript):
315 275
316 def save(self): 276 def save(self):
317 """ 277 """
318 - Save residual map. 278 + Save residual map
319 """ 279 """
320 # Write header 280 # Write header
321 if self._logTerse(): 281 if self._logTerse():
322 - self._log("\n")  
323 - self._log.header1("Save residual map") 282 + self._log('\n')
  283 + self._log.header1('Save residual map')
324 284
325 # Get outmap parameter 285 # Get outmap parameter
326 - outmap = self["outmap"].filename() 286 + outmap = self['outmap'].filename()
327 287
328 # Continue only filename and residual map are valid 288 # Continue only filename and residual map are valid
329 if self._resmap != None: 289 if self._resmap != None:
330 290
331 # Log file name 291 # Log file name
332 if self._logTerse(): 292 if self._logTerse():
333 - self._log(gammalib.parformat("Residual map file")) 293 + self._log(gammalib.parformat('Residual map file'))
334 self._log(outmap.url()) 294 self._log(outmap.url())
335 - self._log("\n") 295 + self._log('\n')
336 296
337 # Save residual map 297 # Save residual map
338 - self._resmap.save(outmap, self._clobber) 298 + self._resmap.save(outmap, self['clobber'].boolean())
339 299
340 # Return 300 # Return
341 return 301 return
342 302
343 - def publish(self, name=""): 303 + def publish(self, name=''):
344 """ 304 """
345 Publish residual map. 305 Publish residual map.
346 306
347 - Kwargs:  
348 - name: Name of residual map. 307 + Parameters
  308 + ----------
  309 + name : str, optional
  310 + Name of residual map
349 """ 311 """
350 # Write header 312 # Write header
351 if self._logTerse(): 313 if self._logTerse():
352 - self._log("\n")  
353 - self._log.header1("Publish residual map") 314 + self._log('\n')
  315 + self._log.header1('Publish residual map')
354 316
355 # Continue only if residual map is valid 317 # Continue only if residual map is valid
356 if self._resmap != None: 318 if self._resmap != None:
357 319
358 # Set default name is user name is empty 320 # Set default name is user name is empty
359 if not name: 321 if not name:
360 - user_name = self.name 322 + user_name = self._name
361 else: 323 else:
362 user_name = name 324 user_name = name
363 325
364 # Log map name 326 # Log map name
365 if self._logTerse(): 327 if self._logTerse():
366 - self._log("Publish residual map \""+user_name+"\".\n") 328 + self._log(gammalib.parformat('Residual map name'))
  329 + self._log(user_name)
  330 + self._log('\n')
367 331
368 # Publish map 332 # Publish map
369 self._resmap.publish(user_name) 333 self._resmap.publish(user_name)
@@ -373,10 +337,10 @@ class csresmap(ctools.cscript): @@ -373,10 +337,10 @@ class csresmap(ctools.cscript):
373 337
374 def models(self, models): 338 def models(self, models):
375 """ 339 """
376 - Set model. 340 + Set model
377 """ 341 """
378 # Copy models 342 # Copy models
379 - self._obs.models(models.clone()) 343 + self._obs.models(models)
380 344
381 # Return 345 # Return
382 return 346 return
test/test_csresmap.py
  View file @63661f5
@@ -76,7 +76,7 @@ class Test(test): @@ -76,7 +76,7 @@ class Test(test):
76 ' emin=0.1 emax=100.0 enumbins=20'+ \ 76 ' emin=0.1 emax=100.0 enumbins=20'+ \
77 ' nxpix=50 nypix=50 binsz=0.02'+ \ 77 ' nxpix=50 nypix=50 binsz=0.02'+ \
78 ' coordsys="CEL" proj="CAR" xref=83.63 yref=22.01'+ \ 78 ' coordsys="CEL" proj="CAR" xref=83.63 yref=22.01'+ \
79 - ' algorithm="SUBDIV"'+ \ 79 + ' algorithm="SUBDIVSQRT"'+ \
80 ' logfile="csresmap_cmd1.log" chatter=1' 80 ' logfile="csresmap_cmd1.log" chatter=1'
81 81
82 # Check if execution of wrong command fails 82 # Check if execution of wrong command fails
@@ -98,7 +98,7 @@ class Test(test): @@ -98,7 +98,7 @@ class Test(test):
98 ' emin=0.1 emax=100.0 enumbins=20'+ \ 98 ' emin=0.1 emax=100.0 enumbins=20'+ \
99 ' nxpix=50 nypix=50 binsz=0.02'+ \ 99 ' nxpix=50 nypix=50 binsz=0.02'+ \
100 ' coordsys="CEL" proj="CAR" xref=83.63 yref=22.01'+ \ 100 ' coordsys="CEL" proj="CAR" xref=83.63 yref=22.01'+ \
101 - ' algorithm="SUBDIV"'+ \ 101 + ' algorithm="SUBDIVSQRT"'+ \
102 ' logfile="csresmap_cmd2.log" chatter=2' 102 ' logfile="csresmap_cmd2.log" chatter=2'
103 103
104 # Check if execution failed 104 # Check if execution failed
@@ -130,7 +130,7 @@ class Test(test): @@ -130,7 +130,7 @@ class Test(test):
130 resmap['proj'] = 'CAR' 130 resmap['proj'] = 'CAR'
131 resmap['xref'] = 83.63 131 resmap['xref'] = 83.63
132 resmap['yref'] = 22.01 132 resmap['yref'] = 22.01
133 - resmap['algorithm'] = 'SUBDIV' 133 + resmap['algorithm'] = 'SUB'
134 resmap['logfile'] = 'csresmap_py1.log' 134 resmap['logfile'] = 'csresmap_py1.log'
135 resmap['chatter'] = 2 135 resmap['chatter'] = 2
136 136
@@ -157,6 +157,7 @@ class Test(test): @@ -157,6 +157,7 @@ class Test(test):
157 resmap['algorithm'] = 'SUBDIV' 157 resmap['algorithm'] = 'SUBDIV'
158 resmap['logfile'] = 'csresmap_py2.log' 158 resmap['logfile'] = 'csresmap_py2.log'
159 resmap['chatter'] = 3 159 resmap['chatter'] = 3
  160 + resmap['publish'] = True
160 161
161 # Run csresmap script 162 # Run csresmap script
162 resmap.execute() 163 resmap.execute()