benchmark_cta_analysis.py
9.93 KB
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
144
145
146
147
148
149
150
151
152
153
154
155
156
157
158
159
160
161
162
163
164
165
166
167
168
169
170
171
172
173
174
175
176
177
178
179
180
181
182
183
184
185
186
187
188
189
190
191
192
193
194
195
196
197
198
199
200
201
202
203
204
205
206
207
208
209
210
211
212
213
214
215
216
217
218
219
220
221
222
223
224
225
226
227
228
229
230
231
232
233
234
235
236
237
238
239
240
241
242
243
244
245
246
247
248
249
250
251
252
253
254
255
256
257
258
259
260
261
262
263
264
265
266
267
268
269
270
271
272
273
274
275
276
277
278
279
280
281
282
283
284
285
286
287
288
289
290
291
292
293
294
295
296
297
298
299
300
301
302
303
304
305
306
307
308
309
310
311
312
313
314
315
316
317
318
319
320
321
322
323
324
325
326
327
328
329
330
331
332
333
334
335
336
337
338
339
340
341
342
343
344
345
346
347
348
349
350
351
352
353
354
355
356
357
358
359
360
361
362
363
#! /usr/bin/env python
# ==========================================================================
# This script peforms a benchmark of the various CTA analysis types
# (unbinned, binned, stacked), and if matplotlib is installed, creates
# a plot of the benchmark results.
#
# Copyright (C) 2014-2015 Jurgen Knodlseder
#
# 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
import ctools
import time
try:
import matplotlib.pyplot as plt
has_matplotlib = True
except:
has_matplotlib = False
# ========================== #
# Unbinned analysis pipeline #
# ========================== #
def unbinned_pipeline(duration):
"""
Unbinned analysis pipeline.
"""
# Set script parameters
model_name = "${CTOOLS}/share/models/crab.xml"
caldb = "prod2"
irf = "South_50h"
ra = 83.63
dec = 22.01
rad_sim = 10.0
tstart = 0.0
tstop = duration
emin = 0.1
emax = 100.0
rad_select = 3.0
# Get start CPU time
tstart = time.clock()
# Simulate events
sim = ctools.ctobssim()
sim["inmodel"].filename(model_name)
sim["caldb"].string(caldb)
sim["irf"].string(irf)
sim["ra"].real(ra)
sim["dec"].real(dec)
sim["rad"].real(rad_sim)
sim["tmin"].real(tstart)
sim["tmax"].real(tstop)
sim["emin"].real(emin)
sim["emax"].real(emax)
sim.run()
# Select events
select = ctools.ctselect(sim.obs())
select["ra"].real(ra)
select["dec"].real(dec)
select["rad"].real(rad_select)
select["tmin"].real(tstart)
select["tmax"].real(tstop)
select["emin"].real(emin)
select["emax"].real(emax)
select.run()
# Get ctlike start CPU time
tctlike = time.clock()
# Perform maximum likelihood fitting
like = ctools.ctlike(select.obs())
like.run()
# Get stop CPU time
tstop = time.clock()
telapsed = tstop - tstart
tctlike = tstop - tctlike
# Return
return telapsed, tctlike
# ======================== #
# Binned analysis pipeline #
# ======================== #
def binned_pipeline(duration):
"""
Binned analysis pipeline.
"""
# Set script parameters
model_name = "${CTOOLS}/share/models/crab.xml"
caldb = "prod2"
irf = "South_50h"
ra = 83.63
dec = 22.01
rad_sim = 10.0
tstart = 0.0
tstop = duration
emin = 0.1
emax = 100.0
enumbins = 20
nxpix = 200
nypix = 200
binsz = 0.02
coordsys = "CEL"
proj = "CAR"
# Get start CPU time
tstart = time.clock()
# Simulate events
sim = ctools.ctobssim()
sim["inmodel"].filename(model_name)
sim["caldb"].string(caldb)
sim["irf"].string(irf)
sim["ra"].real(ra)
sim["dec"].real(dec)
sim["rad"].real(rad_sim)
sim["tmin"].real(tstart)
sim["tmax"].real(tstop)
sim["emin"].real(emin)
sim["emax"].real(emax)
sim.run()
# Bin events into counts map
bin = ctools.ctbin(sim.obs())
bin["ebinalg"].string("LOG")
bin["emin"].real(emin)
bin["emax"].real(emax)
bin["enumbins"].integer(enumbins)
bin["nxpix"].integer(nxpix)
bin["nypix"].integer(nypix)
bin["binsz"].real(binsz)
bin["coordsys"].string(coordsys)
bin["xref"].real(ra)
bin["yref"].real(dec)
bin["proj"].string(proj)
bin.run()
# Get ctlike start CPU time
tctlike = time.clock()
# Perform maximum likelihood fitting
like = ctools.ctlike(bin.obs())
like.run()
# Get stop CPU time
tstop = time.clock()
telapsed = tstop - tstart
tctlike = tstop - tctlike
# Return
return telapsed, tctlike
# ========================= #
# Stacked analysis pipeline #
# ========================= #
def stacked_pipeline(duration):
"""
Cube-style analysis pipeline.
"""
# Set script parameters
model_name = "${CTOOLS}/share/models/crab.xml"
caldb = "prod2"
irf = "South_50h"
ra = 83.63
dec = 22.01
rad_sim = 10.0
tstart = 0.0
tstop = duration
emin = 0.1
emax = 100.0
enumbins = 20
nxpix = 200
nypix = 200
binsz = 0.02
coordsys = "CEL"
proj = "CAR"
# Get start CPU time
tstart = time.clock()
# Simulate events
sim = ctools.ctobssim()
sim["inmodel"].filename(model_name)
sim["caldb"].string(caldb)
sim["irf"].string(irf)
sim["ra"].real(ra)
sim["dec"].real(dec)
sim["rad"].real(rad_sim)
sim["tmin"].real(tstart)
sim["tmax"].real(tstop)
sim["emin"].real(emin)
sim["emax"].real(emax)
sim.run()
# Bin events into counts map
bin = ctools.ctbin(sim.obs())
bin["ebinalg"].string("LOG")
bin["emin"].real(emin)
bin["emax"].real(emax)
bin["enumbins"].integer(enumbins)
bin["nxpix"].integer(nxpix)
bin["nypix"].integer(nypix)
bin["binsz"].real(binsz)
bin["coordsys"].string(coordsys)
bin["proj"].string(proj)
bin["xref"].real(ra)
bin["yref"].real(dec)
bin.run()
# Create exposure cube
expcube = ctools.ctexpcube(sim.obs())
expcube["incube"].filename("NONE")
expcube["caldb"].string(caldb)
expcube["irf"].string(irf)
expcube["ebinalg"].string("LOG")
expcube["emin"].real(emin)
expcube["emax"].real(emax)
expcube["enumbins"].integer(enumbins)
expcube["nxpix"].integer(nxpix)
expcube["nypix"].integer(nypix)
expcube["binsz"].real(binsz)
expcube["coordsys"].string(coordsys)
expcube["proj"].string(proj)
expcube["xref"].real(ra)
expcube["yref"].real(dec)
expcube.run()
# Create PSF cube
psfcube = ctools.ctpsfcube(sim.obs())
psfcube["incube"].filename("NONE")
psfcube["caldb"].string(caldb)
psfcube["irf"].string(irf)
psfcube["ebinalg"].string("LOG")
psfcube["emin"].real(emin)
psfcube["emax"].real(emax)
psfcube["enumbins"].integer(enumbins)
psfcube["nxpix"].integer(10)
psfcube["nypix"].integer(10)
psfcube["binsz"].real(1.0)
psfcube["coordsys"].string(coordsys)
psfcube["proj"].string(proj)
psfcube["xref"].real(ra)
psfcube["yref"].real(dec)
psfcube.run()
# Create background cube
bkgcube = ctools.ctbkgcube(sim.obs())
bkgcube["incube"].filename("NONE")
bkgcube["ebinalg"].string("LOG")
bkgcube["emin"].real(emin)
bkgcube["emax"].real(emax)
bkgcube["enumbins"].integer(enumbins)
bkgcube["nxpix"].integer(10)
bkgcube["nypix"].integer(10)
bkgcube["binsz"].real(1.0)
bkgcube["coordsys"].string(coordsys)
bkgcube["proj"].string(proj)
bkgcube["xref"].real(ra)
bkgcube["yref"].real(dec)
bkgcube.run()
# Attach background model to observation container
bin.obs().models(bkgcube.models())
# Set Exposure and Psf cube for first CTA observation
# (ctbin will create an observation with a single container)
bin.obs()[0].response(expcube.expcube(), psfcube.psfcube(), bkgcube.bkgcube())
# Get ctlike start CPU time
tctlike = time.clock()
# Perform maximum likelihood fitting
like = ctools.ctlike(bin.obs())
like.run()
# Get stop CPU time
tstop = time.clock()
telapsed = tstop - tstart
tctlike = tstop - tctlike
# Return
return telapsed, tctlike
#==========================#
# Main routine entry point #
#==========================#
if __name__ == '__main__':
"""
"""
# Dump header
print("*************************************")
print("* CTA analysis benchmark *")
print("*************************************")
print(" Duration Unbinned Binned Stacked ctlike-UB ctlike-B ctlike-S")
# Initialize arrays
a_duration = []
a_unbinned = []
a_binned = []
a_stacked = []
c_unbinned = []
c_binned = []
c_stacked = []
# Perform benchmarks
duration = 1800.0
for i in range(10):
# Perform analyses
t_unbinned, ct_unbinned = unbinned_pipeline(duration)
t_binned, ct_binned = binned_pipeline(duration)
t_stacked, ct_stacked = stacked_pipeline(duration)
# Collect results
a_duration.append(duration/3600.0)
a_unbinned.append(t_unbinned)
a_binned.append(t_binned)
a_stacked.append(t_stacked)
c_unbinned.append(ct_unbinned)
c_binned.append(ct_binned)
c_stacked.append(ct_stacked)
# Print results
print("%10.0f %10.4f %10.4f %10.4f %10.4f %10.4f %10.4f" % \
(duration, t_unbinned, t_binned, t_stacked, \
ct_unbinned, ct_binned, ct_stacked))
# Increment duration
duration *= 2.0
# Optionally plot results
if has_matplotlib:
plt.figure(1)
plt.title("CTA analysis benchmark")
plt.loglog(a_duration, a_unbinned, 'ro-', label='unbinned')
plt.loglog(a_duration, a_binned, 'bo-', label='binned')
plt.loglog(a_duration, a_stacked, 'go-', label='cube-style')
plt.loglog(a_duration, c_unbinned, 'ro--')
plt.loglog(a_duration, c_binned, 'bo--')
plt.loglog(a_duration, c_stacked, 'go--')
plt.xlabel("Duration (hours)")
plt.ylabel("CPU time (seconds)")
plt.legend(loc="lower right")
plt.show()