-
Notifications
You must be signed in to change notification settings - Fork 10
/
Calibrate.py
executable file
·280 lines (211 loc) · 10.3 KB
/
Calibrate.py
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
#!/usr/bin/env python
# ======================================================================
import pangloss
import sys,getopt,cPickle,numpy
import scipy.stats as stats
# ======================================================================
def Calibrate(argv):
"""
NAME
Calibrate.py
PURPOSE
Transform the results of the lightcone reconstruction process,
Pr(kappah|D), into our target PDF, Pr(kappa|D).
COMMENTS
All PDF input is provided as a list of samples. There are two
modes of operation:
1) The Pr(kappah|C) for an ensemble of calibration lightcones are
compressed into a single number (currently the
median), and then combined with the true kappa values to make
Pr(kappa,kappah|C). This is written out as a 2D sample list.
2) The Pr(kappah|D) for a single observed lightcone is compressed
into a single number (currently the median). This is then used
to take a slice from Pr(kappa,kappah|C) to make Pr(kappa|D,C).
Both 1 and 2 can be carried out in series if desired (Mode=3).
FLAGS
-h Print this message [0]
INPUTS
configfile Plain text file containing Pangloss configuration
OPTIONAL INPUTS
--mode Operating mode 1,2 or 3. See COMMENTS above.
OUTPUTS
stdout Useful information
samples From 1) Pr(kappa,kappah|C) or 2) Pr(kappa|D,C)
EXAMPLE
Calibrate.py example.config
BUGS
AUTHORS
This file is part of the Pangloss project, distributed under the
GPL v2, by Tom Collett (IoA) and Phil Marshall (Oxford).
Please cite: Collett et al 2013, http://arxiv.org/abs/1303.6564
HISTORY
2013-03-21 started Collett & Marshall (Oxford)
"""
# --------------------------------------------------------------------
try:
opts, args = getopt.getopt(argv,"hm:",["help","mode"])
except getopt.GetoptError, err:
print str(err) # will print something like "option -a not recognized"
print Calibrate.__doc__ # will print the big comment above.
return
Mode=3
for o,a in opts:
if o in ("-h", "--help"):
print Calibrate.__doc__
return
elif o in ("-m", "--mode"):
Mode = int(a)
assert Mode < 4 and Mode >0, "unhandled Mode"
else:
assert False, "unhandled option"
# Check for setup file in array args:
if len(args) == 1:
configfile = args[0]
print pangloss.doubledashedline
print pangloss.hello
print pangloss.doubledashedline
print "Calibrate: transforming Pr(kappah|D) to Pr(kappa|D)"
print "Calibrate: taking instructions from",configfile
else:
print Calibrate.__doc__
return
# --------------------------------------------------------------------
# Read in configuration, and extract the ones we need:
experiment = pangloss.Configuration(configfile)
EXP_NAME = experiment.parameters['ExperimentName']
Nc = experiment.parameters['NCalibrationLightcones']
comparator=experiment.parameters['Comparator']
comparatorType=experiment.parameters['ComparatorType']
comparatorWidth=experiment.parameters['ComparatorWidth']
# Figure out which mode is required:
ModeName = experiment.parameters['CalibrateMode']
if ModeName=='Joint': Mode = 1
if ModeName=='Slice': Mode = 2
if ModeName=='JointAndSlice': Mode = 3
CALIB_DIR = experiment.parameters['CalibrationFolder'][0]
jointdistfile= CALIB_DIR+'/'+comparator+'_'+comparatorType+'.pickle'
jointdistasPDFfile= CALIB_DIR+'/'+comparator+'_'+comparatorType+'_asPDF.pickle'
# Final result is PDF for kappa:
x = experiment.parameters['ObservedCatalog'][0]
resultfile = x.split('.')[0]+"_"+EXP_NAME+"_PofKappa.pickle"
# --------------------------------------------------------------------
# Mode 1: generate a joint distribution, eg Pr(kappah,kappa)
# from the calibration dataset:
if Mode==1 or Mode==3:
print pangloss.dashedline
# First find the calibration pdfs for kappa_h:
calpickles = []
for i in range(Nc):
calpickles.append(experiment.getLightconePickleName('simulated',pointing=i))
calresultpickles=[]
if comparator=="Kappah" and comparatorType=="median":
for i in range(Nc):
x = calpickles[i]
pfile = x.split('.')[0].split("_lightcone")[0]+"_"+EXP_NAME+"_KappaHilbert_Kappah_median.pickle"
calresultpickles.append(pfile)
elif comparator=="Kappah" and comparatorType!="median":
for i in range(Nc):
x = calpickles[i]
pfile = x.split('.')[0].split("_lightcone")[0]+"_"+EXP_NAME+"_KappaHilbert_Kappah_"+comparatorType+".pickle"
calresultpickles.append(pfile)
else:
print "Calibrate: Unrecognised comparator "+Comparator
print "Calibrate: If you want to use a comparator other than kappa_h, "
print "Calibrate: you'll need to code it up!"
print "Calibrate: (This should be easy, but you can ask [email protected] for help)."
exit()
# Now calculate comparators:
callist=numpy.empty((Nc,2))
jd=pangloss.PDF(["kappa_ext",comparator+'_'+comparatorType])
for i in range(Nc):
C = calresultpickles[i]
pdf = pangloss.readPickle(C)
if comparator=="Kappah":
if comparatorType=="median":
# Recall that we created a special file for this
# choice of comparator and comparator type, in
# Reconstruct. You could also use the
# comparatortype=="mean" code, swapping mean for median.
callist[i,0]=pdf[0]
callist[i,1]=pdf[1][0]
elif comparatorType=="mean":
callist[i,0] = pdf.truth[0]
callist[i,1] = numpy.mean(pdf.samples)
else:
print "Calibrate: Unrecognised comparatorType "+comparatorType
print "Calibrate: If you want to use a comparatorType other than median "
print "Calibrate: or mean, you'll need to code it up!"
print "Calibrate: (This should be easy, but you can ask [email protected] for help)."
exit()
jd.append(callist[i])
pangloss.writePickle(callist,jointdistfile)
# Also store the joint dist as a pangloss pdf:
pangloss.writePickle(jd,jointdistasPDFfile)
# Plot:
plotfile = jointdistasPDFfile.split('.')[0]+'.png'
jd.plot("Kappah_median","kappa_ext",weight=None,output=plotfile,title="The joint distribution of $\kappa_{\mathrm{ext}}$ and calibrator \n\n (more correlated means a better calibrator!)")
print "Calibrate: calibration joint PDF saved in:"
print "Calibrate: "+jointdistfile
print "Calibrate: and "+jointdistasPDFfile
print "Calibrate: you can view this PDF in "+plotfile
# --------------------------------------------------------------------
# Mode 2: calibrate a real line of sight's Pr(kappah|D) using the
# joint distribution Pr(kappa,<kappah>|D)
if Mode==2 or Mode==3:
print pangloss.dashedline
callibguide = pangloss.readPickle(jointdistfile)
obspickle = experiment.getLightconePickleName('real')
pfile = obspickle.split('.')[0].split("_lightcone")[0]+'_'+EXP_NAME+"_PofKappah.pickle"
pdf=pangloss.readPickle(pfile)
if comparator=="Kappah":
if comparatorType=="median":# note we created a special file for this choice of comparator and comparator type. You could also use the comparatortype=="mean" code swapping mean for median.
RealComparator=numpy.median(pdf.samples)
elif comparatorType=="mean":
RealComparator=numpy.mean(pdf.samples)
else:
print "I don't know that comparatorType. exiting"
exit()
pdf = pangloss.PDF(["kappa_ext","weight"])
#print RealComparator
#print numpy.median(callibguide[:,1]),numpy.std(callibguide[:,1])
dif=(callibguide[:,1]-RealComparator)
weights=dif*0.0
weights[numpy.abs(dif)<comparatorWidth]=1.
weights/=numpy.sum(weights)
samples=callibguide[:,0]
samplesandweights=callibguide.copy()
samplesandweights[:,1]=weights
pdf.samples=(samplesandweights)
plotfile = resultfile.split('.')[0]+".png"
pdf.plot('kappa_ext',weight='weight',output=plotfile)
average = numpy.average(samples, weights=weights)
variance = numpy.dot(weights, (samples-average)**2)/weights.sum()
average,std=average, variance**.5
#if step function weights can calculate 68%CL easily:
included=samples[weights>0]
onesigconfidence=numpy.abs(\
stats.scoreatpercentile(included,84)-
stats.scoreatpercentile(included,16)\
)/2.
pangloss.writePickle(pdf,resultfile)
print "Calibrate: your reconstructed lightcone has been calibrated,"
print "Calibrate: suggesting it has a kappa_ext of",\
"%.3f +\- %.3f"%(average,onesigconfidence)
print "Calibrate: the PDF for kappa_ext has been output to "+resultfile
print "Calibrate: in the form of sample kappa_ext values, and their weights."
print "Calibrate: you can view this PDF in "+plotfile
print
print "Calibrate: To read and process this file, try:"
print
print " import pangloss"
print " pdf = pangloss.readPickle(\"%s\")"%resultfile
print " kappa_samples = pdf.getParameter(\"kappa_ext\")"
print " kappa_weights = pdf.getParameter(\"weight\")"
# --------------------------------------------------------------------
print
print pangloss.doubledashedline
return resultfile,jointdistasPDFfile
# ======================================================================
if __name__ == '__main__':
Calibrate(sys.argv[1:])
# ======================================================================