Mathieu Servillat / ctools

Commit 2bbbeb015632bef16d03364229872314774ef063

Authored by Jürgen Knödlseder
1 parent c9908a82

Add CTA analysis benchmark script

Inline Side-by-side
Showing 1 changed file with 345 additions and 0 deletions
examples/benchmark_cta_analysis.py 0 → 100755
  View file @2bbbeb0
  1 +#! /usr/bin/env python
  2 +# ==========================================================================
  3 +# This script peforms a benchmark of the various CTA analysis styles
  4 +# (unbinned, binned, cube-style).
  5 +#
  6 +# Copyright (C) 2014 Jurgen Knodlseder
  7 +#
  8 +# This program is free software: you can redistribute it and/or modify
  9 +# it under the terms of the GNU General Public License as published by
  10 +# the Free Software Foundation, either version 3 of the License, or
  11 +# (at your option) any later version.
  12 +#
  13 +# This program is distributed in the hope that it will be useful,
  14 +# but WITHOUT ANY WARRANTY; without even the implied warranty of
  15 +# MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
  16 +# GNU General Public License for more details.
  17 +#
  18 +# You should have received a copy of the GNU General Public License
  19 +# along with this program. If not, see <http://www.gnu.org/licenses/>.
  20 +#
  21 +# ==========================================================================
  22 +import gammalib
  23 +import ctools
  24 +import time
  25 +try:
  26 + import matplotlib.pyplot as plt
  27 + has_matplotlib = True
  28 +except:
  29 + has_matplotlib = False
  30 +
  31 +
  32 +# ========================== #
  33 +# Unbinned analysis pipeline #
  34 +# ========================== #
  35 +def unbinned_pipeline(duration):
  36 + """
  37 + Unbinned analysis pipeline.
  38 + """
  39 + # Set script parameters
  40 + model_name = "${CTOOLS}/share/models/crab.xml"
  41 + caldb = "dummy"
  42 + irf = "cta_dummy_irf"
  43 + ra = 83.63
  44 + dec = 22.01
  45 + rad_sim = 10.0
  46 + tstart = 0.0
  47 + tstop = duration
  48 + emin = 0.1
  49 + emax = 100.0
  50 + rad_select = 3.0
  51 +
  52 + # Get start CPU time
  53 + tstart = time.clock()
  54 +
  55 + # Simulate events
  56 + sim = ctools.ctobssim()
  57 + sim["infile"].filename(model_name)
  58 + sim["caldb"].string(caldb)
  59 + sim["irf"].string(irf)
  60 + sim["ra"].real(ra)
  61 + sim["dec"].real(dec)
  62 + sim["rad"].real(rad_sim)
  63 + sim["tmin"].real(tstart)
  64 + sim["tmax"].real(tstop)
  65 + sim["emin"].real(emin)
  66 + sim["emax"].real(emax)
  67 + sim.run()
  68 +
  69 + # Select events
  70 + select = ctools.ctselect(sim.obs())
  71 + select["ra"].real(ra)
  72 + select["dec"].real(dec)
  73 + select["rad"].real(rad_select)
  74 + select["tmin"].real(tstart)
  75 + select["tmax"].real(tstop)
  76 + select["emin"].real(emin)
  77 + select["emax"].real(emax)
  78 + select.run()
  79 +
  80 + # Get ctlike start CPU time
  81 + tctlike = time.clock()
  82 +
  83 + # Perform maximum likelihood fitting
  84 + like = ctools.ctlike(select.obs())
  85 + like.run()
  86 +
  87 + # Get stop CPU time
  88 + tstop = time.clock()
  89 + telapsed = tstop - tstart
  90 + tctlike = tstop - tctlike
  91 +
  92 + # Return
  93 + return telapsed, tctlike
  94 +
  95 +
  96 +# ======================== #
  97 +# Binned analysis pipeline #
  98 +# ======================== #
  99 +def binned_pipeline(duration):
  100 + """
  101 + Binned analysis pipeline.
  102 + """
  103 + # Set script parameters
  104 + model_name = "${CTOOLS}/share/models/crab.xml"
  105 + caldb = "dummy"
  106 + irf = "cta_dummy_irf"
  107 + ra = 83.63
  108 + dec = 22.01
  109 + rad_sim = 10.0
  110 + tstart = 0.0
  111 + tstop = duration
  112 + emin = 0.1
  113 + emax = 100.0
  114 + enumbins = 20
  115 + nxpix = 200
  116 + nypix = 200
  117 + binsz = 0.02
  118 + coordsys = "CEL"
  119 + proj = "CAR"
  120 +
  121 + # Get start CPU time
  122 + tstart = time.clock()
  123 +
  124 + # Simulate events
  125 + sim = ctools.ctobssim()
  126 + sim["infile"].filename(model_name)
  127 + sim["caldb"].string(caldb)
  128 + sim["irf"].string(irf)
  129 + sim["ra"].real(ra)
  130 + sim["dec"].real(dec)
  131 + sim["rad"].real(rad_sim)
  132 + sim["tmin"].real(tstart)
  133 + sim["tmax"].real(tstop)
  134 + sim["emin"].real(emin)
  135 + sim["emax"].real(emax)
  136 + sim.run()
  137 +
  138 + # Bin events into counts map
  139 + bin = ctools.ctbin(sim.obs())
  140 + bin["ebinalg"].string("LOG")
  141 + bin["emin"].real(emin)
  142 + bin["emax"].real(emax)
  143 + bin["enumbins"].integer(enumbins)
  144 + bin["nxpix"].integer(nxpix)
  145 + bin["nypix"].integer(nypix)
  146 + bin["binsz"].real(binsz)
  147 + bin["coordsys"].string(coordsys)
  148 + bin["xref"].real(ra)
  149 + bin["yref"].real(dec)
  150 + bin["proj"].string(proj)
  151 + bin.run()
  152 +
  153 + # Get ctlike start CPU time
  154 + tctlike = time.clock()
  155 +
  156 + # Perform maximum likelihood fitting
  157 + like = ctools.ctlike(bin.obs())
  158 + like.run()
  159 +
  160 + # Get stop CPU time
  161 + tstop = time.clock()
  162 + telapsed = tstop - tstart
  163 + tctlike = tstop - tctlike
  164 +
  165 + # Return
  166 + return telapsed, tctlike
  167 +
  168 +
  169 +# ============================ #
  170 +# Cube-style analysis pipeline #
  171 +# ============================ #
  172 +def cube_pipeline(duration):
  173 + """
  174 + Cube-style analysis pipeline.
  175 + """
  176 + # Set script parameters
  177 + model_name = "${CTOOLS}/share/models/crab.xml"
  178 + caldb = "dummy"
  179 + irf = "cta_dummy_irf"
  180 + ra = 83.63
  181 + dec = 22.01
  182 + rad_sim = 10.0
  183 + tstart = 0.0
  184 + tstop = duration
  185 + emin = 0.1
  186 + emax = 100.0
  187 + enumbins = 20
  188 + nxpix = 200
  189 + nypix = 200
  190 + binsz = 0.02
  191 + coordsys = "CEL"
  192 + proj = "CAR"
  193 +
  194 + # Get start CPU time
  195 + tstart = time.clock()
  196 +
  197 + # Simulate events
  198 + sim = ctools.ctobssim()
  199 + sim["infile"].filename(model_name)
  200 + sim["caldb"].string(caldb)
  201 + sim["irf"].string(irf)
  202 + sim["ra"].real(ra)
  203 + sim["dec"].real(dec)
  204 + sim["rad"].real(rad_sim)
  205 + sim["tmin"].real(tstart)
  206 + sim["tmax"].real(tstop)
  207 + sim["emin"].real(emin)
  208 + sim["emax"].real(emax)
  209 + sim.run()
  210 +
  211 + # Bin events into counts map
  212 + bin = ctools.ctbin(sim.obs())
  213 + bin["ebinalg"].string("LOG")
  214 + bin["emin"].real(emin)
  215 + bin["emax"].real(emax)
  216 + bin["enumbins"].integer(enumbins)
  217 + bin["nxpix"].integer(nxpix)
  218 + bin["nypix"].integer(nypix)
  219 + bin["binsz"].real(binsz)
  220 + bin["coordsys"].string(coordsys)
  221 + bin["proj"].string(proj)
  222 + bin["xref"].real(ra)
  223 + bin["yref"].real(dec)
  224 + bin.run()
  225 +
  226 + # Create exposure cube
  227 + expcube = ctools.ctexpcube(bin.obs())
  228 + expcube["cntmap"].filename("NONE")
  229 + expcube["caldb"].string(caldb)
  230 + expcube["irf"].string(irf)
  231 + expcube["ebinalg"].string("LOG")
  232 + expcube["emin"].real(emin)
  233 + expcube["emax"].real(emax)
  234 + expcube["enumbins"].integer(enumbins)
  235 + expcube["nxpix"].integer(nxpix)
  236 + expcube["nypix"].integer(nypix)
  237 + expcube["binsz"].real(binsz)
  238 + expcube["coordsys"].string(coordsys)
  239 + expcube["proj"].string(proj)
  240 + expcube["xref"].real(ra)
  241 + expcube["yref"].real(dec)
  242 + expcube.run()
  243 +
  244 + # Create PSF cube
  245 + psfcube = ctools.ctpsfcube(bin.obs())
  246 + psfcube["cntmap"].filename("NONE")
  247 + psfcube["caldb"].string(caldb)
  248 + psfcube["irf"].string(irf)
  249 + psfcube["ebinalg"].string("LOG")
  250 + psfcube["emin"].real(emin)
  251 + psfcube["emax"].real(emax)
  252 + psfcube["enumbins"].integer(enumbins)
  253 + psfcube["nxpix"].integer(10)
  254 + psfcube["nypix"].integer(10)
  255 + psfcube["binsz"].real(1.0)
  256 + psfcube["coordsys"].string(coordsys)
  257 + psfcube["proj"].string(proj)
  258 + psfcube["xref"].real(ra)
  259 + psfcube["yref"].real(dec)
  260 + psfcube.run()
  261 +
  262 + # Set exposure and PSF cube for first CTA observation
  263 + obs = bin.obs()
  264 + cta = gammalib.GCTAObservation(obs[0])
  265 + cta.response(expcube.expcube(), psfcube.psfcube())
  266 + obs[0] = cta
  267 +
  268 + # Get ctlike start CPU time
  269 + tctlike = time.clock()
  270 +
  271 + # Perform maximum likelihood fitting
  272 + like = ctools.ctlike(obs)
  273 + like.run()
  274 +
  275 + # Get stop CPU time
  276 + tstop = time.clock()
  277 + telapsed = tstop - tstart
  278 + tctlike = tstop - tctlike
  279 +
  280 + # Return
  281 + return telapsed, tctlike
  282 +
  283 +
  284 +#==========================#
  285 +# Main routine entry point #
  286 +#==========================#
  287 +if __name__ == '__main__':
  288 + """
  289 + """
  290 + # Dump header
  291 + print("*************************************")
  292 + print("* CTA analysis benchmark *")
  293 + print("*************************************")
  294 + print(" Duration Unbinned Binned Cube-style ctlike-UB ctlike-B ctlike-CS")
  295 +
  296 + # Initialize arrays
  297 + a_duration = []
  298 + a_unbinned = []
  299 + a_binned = []
  300 + a_cube = []
  301 + c_unbinned = []
  302 + c_binned = []
  303 + c_cube = []
  304 +
  305 + # Perform benchmarks
  306 + duration = 1800.0
  307 + for i in range(10):
  308 +
  309 + # Perform analyses
  310 + t_unbinned, ct_unbinned = unbinned_pipeline(duration)
  311 + t_binned, ct_binned = binned_pipeline(duration)
  312 + t_cube, ct_cube = cube_pipeline(duration)
  313 +
  314 + # Collect results
  315 + a_duration.append(duration/3600.0)
  316 + a_unbinned.append(t_unbinned)
  317 + a_binned.append(t_binned)
  318 + a_cube.append(t_cube)
  319 + c_unbinned.append(ct_unbinned)
  320 + c_binned.append(ct_binned)
  321 + c_cube.append(ct_cube)
  322 +
  323 + # Print results
  324 + print("%10.0f %10.4f %10.4f %10.4f %10.4f %10.4f %10.4f" % \
  325 + (duration, t_unbinned, t_binned, t_cube, \
  326 + ct_unbinned, ct_binned, ct_cube))
  327 +
  328 + # Increment duration
  329 + duration *= 2.0
  330 +
  331 + # Optionally plot results
  332 + if has_matplotlib:
  333 + plt.figure(1)
  334 + plt.title("CTA analysis benchmark")
  335 + plt.loglog(a_duration, a_unbinned, 'ro-', label='unbinned')
  336 + plt.loglog(a_duration, a_binned, 'bo-', label='binned')
  337 + plt.loglog(a_duration, a_cube, 'go-', label='cube-style')
  338 + plt.loglog(a_duration, c_unbinned, 'ro--')
  339 + plt.loglog(a_duration, c_binned, 'bo--')
  340 + plt.loglog(a_duration, c_cube, 'go--')
  341 + plt.xlabel("Duration (hours)")
  342 + plt.ylabel("CPU time (seconds)")
  343 + plt.legend(loc="lower right")
  344 + plt.show()
  345 +
0 \ No newline at end of file 346 \ No newline at end of file