-
Notifications
You must be signed in to change notification settings - Fork 1
/
Copy pathquadoctcorr.py
298 lines (239 loc) · 10.1 KB
/
quadoctcorr.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
281
282
283
284
285
286
287
288
289
290
291
292
293
294
295
296
#! /usr/bin/env python
"""
Name:
quadoctcorr (quadrupole-octopole and correlation)
Purpose:
extract l=2,3 a_lm.s from SMICA to compare C(theta) against simulations
Uses:
healpy
legprodint
get_crosspower
sim_stats
ispice
Inputs:
Outputs:
Modification History:
Largely copied from sim_stats; Z Knight, 2016.09.13
Added [:lmax+1] to getCovar calls; ZK, 2016.09.14
Switched to useLensing=1 in loadCls; ZK, 2016.10.07
"""
import numpy as np
import matplotlib.pyplot as plt
import healpy as hp
import time # for measuring duration
from scipy.interpolate import interp1d
import get_crosspower as gcp # for loadCls
from numpy.polynomial.legendre import legval # for C_l -> C(theta) conversion
from numpy.polynomial.legendre import legfit # for C(theta) -> C_l conversion
from ispice import ispice
from legprodint import getJmn
import sim_stats as sims # for getSMICA and getCovar
def mf():
"""
Purpose:
Procedure:
Inputs:
Returns:
"""
pass
################################################################################
# testing code
def test(useCLASS=1,useLensing=1,classCamb=1,nSims=1000,lmax=3,lmin=2,
newSMICA=True,newDeg=False):
"""
code for testing the other functions in this module
Inputs:
useCLASS: set to 1 to use CLASS, 0 to use CAMB
CLASS Cl has early/late split at z=50
CAMB Cl has ISWin/out split: ISWin: 0.4<z<0.75, ISWout: the rest
Note: CAMB results include primary in ISWin and ISWout (not as intended)
default: 1
useLensing: set to 1 to use lensed Cl, 0 for non-lensed
default: 1
classCamb: if 1: use the CAMB format of CLASS output, if 0: use CLASS format
Note: parameter not used if useCLASS = 0
default: 1
nSims: the number of simulations to do for ensemble
default: 1000
lmax: the highest l to include in Legendre transforms
default: 3
lmin: the lowest l to include in S_{1/2} = CIC calculations
default: 2
newSMICA: set to True to recalculate SMICA results
default: True
newDeg: set to True to recalculate map and mask degredations
default: False
"""
##############################################################################
# load theoretical power spectra
# load data
ell,fullCl,primCl,lateCl,crossCl = gcp.loadCls(useCLASS=useCLASS,useLensing=useLensing,classCamb=classCamb)
# fill beginning with zeros
startEll = ell[0]
ell = np.append(np.arange(startEll),ell)
fullCl = np.append(np.zeros(startEll),fullCl)
primCl = np.append(np.zeros(startEll),primCl)
lateCl = np.append(np.zeros(startEll),lateCl)
crossCl = np.append(np.zeros(startEll),crossCl)
# suppress C_2 to see what happens in enesmble
suppressC2 = False
suppFactor = 0.23 # from Tegmark et. al. 2003, figure 13 (WMAP)
if suppressC2:
fullCl[2] *= suppFactor
primCl[2] *= suppFactor
lateCl[2] *= suppFactor
crossCl[2] *= suppFactor
conv = ell*(ell+1)/(2*np.pi)
#print ell,conv #ell[0]=2.0
# apply beam and pixel window functions to power spectra
# note: to ignore the non-constant pixel shape, W(l) must be > B(l)
# however, this is not true for NSIDE=128 and gauss_beam(5')
# Here I ignore this anyway and proceed
myNSIDE = 128 # must be same NSIDE as in sims.getSMICA function
Wpix = hp.pixwin(myNSIDE)
Bsmica = hp.gauss_beam(5./60*np.pi/180) # 5 arcmin
WlMax = Wpix.size
if WlMax < lmax:
print 'die screaming!!!'
return 0
fullCl = fullCl[:WlMax]*(Wpix*Bsmica)**2
primCl = primCl[:WlMax]*(Wpix*Bsmica)**2
lateCl = lateCl[:WlMax]*(Wpix*Bsmica)**2
crossCl = crossCl[:WlMax]*(Wpix*Bsmica)**2
# note: i tried sims without this scaling, and results seemed the same at a glance
# extract the part I want
myL = ell[:lmax]
myCl = fullCl[:lmax]
##############################################################################
# load SMICA data and filter out all but low-l a_lm.s
theta_i = 0.0 #degrees
theta_f = 180.0 #degrees
nSteps = 1800
# default: lmax=3,lmin=2
#newSMICA = True # so I don't use lmax=100 from previous calc.
thetaArray2sp, C_SMICAsp, C_SMICAmaskedsp, S_SMICAnomasksp, S_SMICAmaskedsp = \
sims.getSMICA(theta_i=theta_i,theta_f=theta_f,nSteps=nSteps,lmax=lmax,lmin=lmin,
newSMICA=newSMICA,newDeg=newDeg,useSPICE=True)
##############################################################################
# create ensemble of realizations and gather statistics
covEnsembleFull = np.zeros([nSims,nSteps+1]) # for maskless
covEnsembleCut = np.zeros([nSims,nSteps+1]) # for masked
sEnsembleFull = np.zeros(nSims)
sEnsembleCut = np.zeros(nSims)
covTheta = np.array([])
# get Jmn matrix for harmonic space S_{1/2} calc.
myJmn = getJmn(lmax=lmax)
doTime = True # to time the run and print output
startTime = time.time()
for nSim in range(nSims):
print 'starting sim ',nSim+1, ' of ',nSims
alm_prim,alm_late = hp.synalm((primCl,lateCl,crossCl),lmax=lmax,new=True)
# calculate C(theta) of simulation
Clsim_prim = hp.alm2cl(alm_prim)
Clsim_late = hp.alm2cl(alm_late)
Clsim_cros = hp.alm2cl(alm_prim,alm_late)
Clsim_full = Clsim_prim + 2*Clsim_cros + Clsim_late
# use Cl_sim_full to omit prim/late distinction for now
#Clsim_full_sum += Clsim_full
# first without mask
# note: getCovar uses linspace in x for thetaArray
thetaArray,cArray = sims.getCovar(ell[:lmax+1],Clsim_full[:lmax+1],theta_i=theta_i,
theta_f=theta_f,nSteps=nSteps,lmin=lmin)
covEnsembleFull[nSim] = cArray
covTheta = thetaArray
# S_{1/2}
sEnsembleFull[nSim] = np.dot(Clsim_full[lmin:],np.dot(myJmn[lmin:,lmin:],Clsim_full[lmin:]))
# now with a mask
# should have default RING ordering
# pixel window and beam already accounted for in true Cls
#mapSim = hp.alm2map(alm_prim+alm_late,myNSIDE,lmax=lmax,pixwin=True,sigma=5./60*np.pi/180)
mapSim = hp.alm2map(alm_prim+alm_late,myNSIDE,lmax=lmax)
# super lame that spice needs to read/write from disk, but here goes...
mapTempFile = 'tempMap.fits'
ClTempFile = 'tempCl.fits'
maskDegFile = 'maskMapDeg.fits' # this should have been created by sims.getSMICA
hp.write_map(mapTempFile,mapSim)
ispice(mapTempFile,ClTempFile,maskfile1=maskDegFile,subav="YES",subdipole="YES")
Cl_masked = hp.read_cl(ClTempFile)
ell2 = np.arange(Cl_masked.shape[0])
# note: getCovar uses linspace in x for thetaArray
thetaArray,cArray2 = sims.getCovar(ell2[:lmax+1],Cl_masked[:lmax+1],theta_i=theta_i,
theta_f=theta_f,nSteps=nSteps,lmin=lmin)
covEnsembleCut[nSim] = cArray2
# S_{1/2}
sEnsembleCut[nSim] = np.dot(Cl_masked[lmin:lmax+1],np.dot(myJmn[lmin:,lmin:],
Cl_masked[lmin:lmax+1]))
doPlot = False#True
if doPlot:
plt.plot(thetaArray,cArray)
plt.xlabel('theta (degrees)')
plt.ylabel('C(theta)')
plt.title('covariance of simulated CMB')
plt.show()
if doTime: print 'time elapsed: ',int((time.time()-startTime)/60.),' minutes'
avgEnsembleFull = np.average(covEnsembleFull, axis = 0)
stdEnsembleFull = np.std(covEnsembleFull, axis = 0)
# do I need a better way to describe confidence interval?
avgEnsembleCut = np.average(covEnsembleCut, axis = 0)
stdEnsembleCut = np.std(covEnsembleCut, axis = 0)
#Clsim_full_avg = Clsim_full_sum / nSims
##############################################################################
# plot/print results
doPlot = True
if doPlot:
print 'plotting correlation functions... '
# first the whole sky statistics
plt.plot(thetaArray,avgEnsembleFull,label='sim. ensemble average (no mask)')
plt.fill_between(thetaArray,avgEnsembleFull+stdEnsembleFull,
avgEnsembleFull-stdEnsembleFull,alpha=0.25,
label='simulation 1sigma envelope')
#plt.plot(thetaArray2,C_SMICA,label='SMICA R2 (inpainted,anafast)')
plt.plot(thetaArray2sp,C_SMICAsp,label='SMICA R2 (inpainted,spice)')
plt.xlabel('theta (degrees)')
plt.ylabel('C(theta)')
plt.title('whole sky covariance of '+str(nSims)+' simulated CMBs, lmax='+str(lmax))
plt.ylim([-500,1000])
plt.plot([0,180],[0,0]) #horizontal line
plt.legend()
plt.show()
# now the cut sky
plt.plot(thetaArray,avgEnsembleCut,label='sim. ensemble average (masked)')
plt.fill_between(thetaArray,avgEnsembleCut+stdEnsembleCut,
avgEnsembleCut-stdEnsembleCut,alpha=0.25,
label='simulation 1sigma envelope')
#plt.plot(thetaArray2,C_SMICAmasked,label='SMICA R2 (masked ,anafast)')
plt.plot(thetaArray2sp,C_SMICAmaskedsp,label='SMICA R2 (masked ,spice)')
plt.xlabel('theta (degrees)')
plt.ylabel('C(theta)')
plt.title('cut sky covariance of '+str(nSims)+' simulated CMBs, lmax='+str(lmax))
plt.ylim([-500,1000])
plt.plot([0,180],[0,0]) #horizontal line
plt.legend()
plt.show()
print 'plotting S_{1/2} distributions... '
myBins = np.logspace(2,7,100)
plt.axvline(x=S_SMICAnomasksp,color='b',linewidth=3,label='SMICA inpainted')
plt.axvline(x=S_SMICAmaskedsp,color='g',linewidth=3,label='SMICA masked')
plt.hist(sEnsembleFull,bins=myBins,histtype='step',label='full sky')
plt.hist(sEnsembleCut, bins=myBins,histtype='step',label='cut sky')
plt.gca().set_xscale("log")
plt.legend()
plt.xlabel('S_{1/2} (microK^4)')
plt.ylabel('Counts')
plt.title('S_{1/2} of '+str(nSims)+' simulated CMBs')
plt.show()
# S_{1/2} output
print ''
print 'using CIC method: '
#print 'S_{1/2}(anafast): SMICA, no mask: ',S_SMICAnomask,', masked: ',S_SMICAmasked
print 'S_{1/2}(spice): SMICA, no mask: ',S_SMICAnomasksp,', masked: ',S_SMICAmaskedsp
print ''
#print 'using CCdx method: '
#print 'S_{1/2}(anafast): SMICA, no mask: ',SSnm2,', masked: ',SSmd2
#print 'S_{1/2}(spice): SMICA, no mask: ',SSnm2sp,', masked: ',SSmd2sp
#print ''
##############################################################################
# step 3
print 'step 3: profit'
if __name__=='__main__':
test()