Commit 8e15f20ece051f1004aeebb06d465a4f2384c430
1 parent
78d60905
removed (replaced by cspull script)
Showing
1 changed file
with
0 additions
and
141 deletions
test/test_coverage.py deleted
| 1 | -#! /usr/bin/env python | |
| 2 | -# ========================================================================== | |
| 3 | -# This script computes the pull distribution for models. | |
| 4 | -# | |
| 5 | -# Copyright (C) 2011 Christoph Deil | |
| 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 | -from ctatools import * | |
| 22 | -from gammalib import * | |
| 23 | -from math import * | |
| 24 | -import os | |
| 25 | -import glob | |
| 26 | -import sys | |
| 27 | -import csv | |
| 28 | - | |
| 29 | -# Parameters | |
| 30 | -events_name = "events.fits" | |
| 31 | -cntmap_name = "cntmap.fits" | |
| 32 | -result_name = "results.xml" | |
| 33 | -caldb = "irf" | |
| 34 | -irf = "kb_E_50h_v3" | |
| 35 | -ra = 83.63 | |
| 36 | -dec = 22.01 | |
| 37 | -rad_sim = 10.0 | |
| 38 | -tstart = 0.0 | |
| 39 | -tstop = 1e4 | |
| 40 | -emin = 1 | |
| 41 | -emax = 10 | |
| 42 | -enumbins = 5 | |
| 43 | -nxpix = 200 | |
| 44 | -nypix = 200 | |
| 45 | -binsz = 0.02 | |
| 46 | -coordsys = "CEL" | |
| 47 | -proj = "CAR" | |
| 48 | - | |
| 49 | -def run_sim(seed, model_name): | |
| 50 | - sim = ctobssim() | |
| 51 | - sim.logFileOpen() # We need this to explicitely open the log file in Python mode | |
| 52 | - sim["infile"].filename(model_name) | |
| 53 | - sim["outfile"].filename(events_name) | |
| 54 | - sim["caldb"].string(caldb) | |
| 55 | - sim["irf"].string(irf) | |
| 56 | - sim["seed"].integer(seed) | |
| 57 | - sim["ra"].real(ra) | |
| 58 | - sim["dec"].real(dec) | |
| 59 | - sim["rad"].real(rad_sim) | |
| 60 | - sim["tmin"].real(tstart) | |
| 61 | - sim["tmax"].real(tstop) | |
| 62 | - sim["emin"].real(emin) | |
| 63 | - sim["emax"].real(emax) | |
| 64 | - sim.execute() | |
| 65 | - return sim | |
| 66 | - | |
| 67 | -def run_bin(): | |
| 68 | - bin = ctbin() | |
| 69 | - bin.logFileOpen() # We need this to explicitely open the log file in Python mode | |
| 70 | - bin["evfile"].filename(events_name) | |
| 71 | - bin["outfile"].filename(cntmap_name) | |
| 72 | - bin["emin"].real(emin) | |
| 73 | - bin["emax"].real(emax) | |
| 74 | - bin["enumbins"].integer(enumbins) | |
| 75 | - bin["nxpix"].integer(nxpix) | |
| 76 | - bin["nypix"].integer(nypix) | |
| 77 | - bin["binsz"].real(binsz) | |
| 78 | - bin["coordsys"].string(coordsys) | |
| 79 | - bin["xref"].real(ra) | |
| 80 | - bin["yref"].real(dec) | |
| 81 | - bin["proj"].string(proj) | |
| 82 | - bin.execute() | |
| 83 | - return bin | |
| 84 | - | |
| 85 | -def run_like(model_name): | |
| 86 | - like = ctlike() | |
| 87 | - like.logFileOpen() # We need this to explicitely open the log file in Python mode | |
| 88 | - like["cntmap"].filename(cntmap_name) | |
| 89 | - like["srcmdl"].filename(model_name) | |
| 90 | - like["outmdl"].filename(result_name) | |
| 91 | - like["method"].string("BINNED") | |
| 92 | - like["caldb"].string(caldb) | |
| 93 | - like["irf"].string(irf) | |
| 94 | - like.execute() | |
| 95 | - return like | |
| 96 | - | |
| 97 | -def test_coverage(n_trials=3, model_name="data/crab.xml"): | |
| 98 | - """Simulate and fit ntrials time and look at the | |
| 99 | - pull distribution for all free parameters. | |
| 100 | - It should be Gaussian with mean 0 and rms 1""" | |
| 101 | - | |
| 102 | - results = [] | |
| 103 | - colnames = [] # to preserve the order! | |
| 104 | - for seed in range(n_trials): | |
| 105 | - print 'seed:', seed | |
| 106 | - sim = run_sim(seed, model_name) | |
| 107 | - run_bin() | |
| 108 | - like = run_like(model_name) | |
| 109 | - models = like.obs().models() | |
| 110 | - result = {} | |
| 111 | - colnames = [] | |
| 112 | - for par in range(models.npars()): | |
| 113 | - p = models.par(par) | |
| 114 | - if p.isfree(): | |
| 115 | - name = p.name() | |
| 116 | - colnames.append(name) | |
| 117 | - fitted_value = p.value() | |
| 118 | - real_value = sim.obs().models().par(par).value() | |
| 119 | - error = p.error() | |
| 120 | - pull = (fitted_value - real_value) / error | |
| 121 | - result[name] = pull | |
| 122 | - results.append(result) | |
| 123 | - return results, colnames | |
| 124 | - | |
| 125 | -if __name__ == '__main__': | |
| 126 | - usage = "test_coverage xml_model_file n_trials results_file" | |
| 127 | - if len(sys.argv) < 4: | |
| 128 | - print usage | |
| 129 | - sys.exit() | |
| 130 | - | |
| 131 | - xml_model_file = sys.argv[1] | |
| 132 | - n_trials = int(sys.argv[2]) | |
| 133 | - results_file = sys.argv[3] | |
| 134 | - | |
| 135 | - results, colnames = test_coverage(n_trials, xml_model_file) | |
| 136 | - | |
| 137 | - f = open(results_file, 'w') | |
| 138 | - writer = csv.DictWriter(f, colnames) | |
| 139 | - writer.writerow(dict((_,_) for _ in colnames)) | |
| 140 | - writer.writerows(results) | |
| 141 | - f.close() | |
| 142 | 0 | \ No newline at end of file |