1:
HWK3
November 18, 2019
1 Planck 1d circular slices data project
Note: Im using weave to make a notebook from a script file.
using Weave
convertdocnotebook2.jl, notebook2.md for markdown format convertdocnotebook2.jl, notebook2.ipynb for jupyter notebook format
Or directly from shell
julia e using Weave; convertdocnotebook2.jl, notebook2.ipynb
julia e using Weave; convertdocnotebook2.jl, notebook2.md
1.1 Start by defining a module to hold some constants and set the path to the data directory.
module LocalMethods
using GSL
using LinearAlgebra using Statistics import JLD2
export covn1dotn2, covangle
const datadir UsersethananderesDropboxCoursesSTA250fall2019Datasets Dataset2Planck1dtransectdata
const cls JLD2.FileIO.loaddatadircls.jld2, cls
const wonPl cls:cllenscalar:,1 . 2 . cls:ell . 1 . 4 .
cls:factoronclcmb
const lmax maximumcls:ell
function covcos cos ,lmax::Int
assert lmax 2 lmax too small, must be greater than 1 dotsflegendrePlarraylmax, cos 3:end, wonPl3:lmax1
end
covn1dotn2n1dotn2 covcos n1dotn2,lmax
1
covangleangle covcos cosangle,lmax
function binmeanfk; bin0 fkrtn copyfk
if bin 1
Nm1 lengthfk
subNm1 bin Nm1bin
fmk reshapefk1:subNm1, bin, Nm1bin fmk . meanfmk, dims1
fkrtn1:subNm1 . vecfmk fkrtnsubNm11:end . meanfksubNm11:end
end
return fkrtn
end
end module
1: Main.LocalMethods 1.2 Load modules
2:
using PyPlot
using PyCall
using LinearAlgebra using Statistics using Random import JLD2
using FFTW
using SparseArrays using ProgressMeter using .LocalMethods
const propcycle plt.rcParamsaxes.propcycle const colors propcycle.bykeycolor const mplot3d pyimportmpltoolkits.mplot3d
PyObject module mpltoolkits.mplot3d from
UsersethananderesSoftwareanaconda3libpython3.7site
packagesmpltoolkitsmplot3dinit.py
2:
1.3 Load the data
The data loaded below contains 10 slices of the CMB at 10 different polar angles. 5 slices from the north pole, 5 slices from the south pole. The polar angles from which these slices are taken are given in variables northp and southp. The actual CMB measurements at each slice is stored
2
in variables cmbnorthp and cmbsouthp which are simply vectors the indices of which map to azmuthal coordinates stored in variable ref.
These slices are all taken at polar angles which have the same azmuthal coordinates, i.e. each CMB slice is the same length denoted N below the entries of which have asmuthal coordinates given in
ref. 3:
loadit jld2file, varname JLD2.FileIO.loadLocalMethods. datadirjld2file, varname
northp loaditdataslices.jld2, northp
southp loaditdataslices.jld2, southp cmbnorthp loaditdataslices.jld2, cmbnorthp cmbsouthp loaditdataslices.jld2, cmbsouthp
ref loaditdataslices.jld2, ref N length ref
Nside 2048
3: 2048
For the homework assignemnt, only analyze the first and fourth CMB slice from the south pole.
Here is the code for extracting them.
4:
Islice1 1
Islice2 4
1 southpIslice1
2 southpIslice2 cmb1 cmbsouthpIslice1 cmb2 cmbsouthpIslice2
show 1, 2
show summarycmb1 show summarycmb2 show summary ref;
1, 2 1.990924632438072, 2.2330667085254525 summarycmb1 8192element ArrayFloat64,1 summarycmb2 8192element ArrayFloat64,1 summary ref 8192element ArrayFloat64,1
ofsetcmb9000.vcatnorthp,southp. .0.5reverse figurefigsize10,5
for i 1:length northp
plot ref, cmbnorthpi . ofsetcmbi, k, alpha 0.1
end
for i 1:length southp
if i in Islice1, Islice2
plot ref, cmbsouthpi . ofsetcmbilength southp
else
5:
3
alpha 0.1 end
end
plot ref, cmbsouthpi . ofsetcmbilength southp, k,
Note:ThecmbisinunitsmK MicroKelvin1Kelvin1e6MicroKelvin
1.4 Digression: Take a peak at the full sky data
The full sky data file is also uploaded to Canvas file datadatafull.jld2 but you dont need it unless you want to look at it. In which case you will need healpy to make the following few plots. install healpy at the command line with something like pip install user healpy
const hp pyimporthealpy
const ScS pyimportscipy.spatial const ScI pyimportscipy.interpolate
cmbfull loaditdatafull.jld2, cmbfull full loaditdatafull.jld2, full full loaditdatafull.jld2, full
hp.visufunc.mollviewcmbfull, xsize2048, titlefull cmb map
hp.visufunc.mollviewcmbfull, min250, max500, xsizeNside, titlefull cmb map with restricted color scale
6:
4
To get a visual where the CMB slices used in the homework come from, they are take at polar angles from within these strips.
5
custommask northp1 . full . northpend . southp1 . full . southpend
hp.visufunc.mollviewcustommask . cmbfull, min250, max500, xsizeNside
7:
1.5 1st column of covariance covariance and crosscovariance matrix
The isotropic nature of the CMB means we can compute the covariance between the pixel obser vations as a function of the inner product of the three dimensional sphere coordinates. Here is a function that eats polar coordaintes and returns the tuple of three dim spherical coordinates.
function n ,
nx sin..cos.
ny sin..sin. nz cos.
return nx, ny, nz
end
show n 1, ref1; show n 2, ref1;
8:
n 1, ref1 0.9130366148294711, 0.0, 0.4078776041666667 n 2, ref1 0.7885981873346127, 0.0, 0.6149088541666667
Now we can compute the inner product between a fixed pixel and all the rest on the same slice as follows
6
we use Ref here since n 1, ref1 is a single 3tuple that we want broadcasted to the vector of 3tuples n. 1, ref
dot.n.1, ref,Refn1, ref1
9:
9:
8192element ArrayFloat64,1:
1.0
0.9999997547967509
0.9999990191871476
0.9999977931716232
0.9999960767508986
0.9999938699259837
0.9999911726981767
0.9999879850690642
0.9999843070405215
0.9999801386147124
0.9999754797940889
0.9999703305813916
0.9999646909796499
0.9999646909796499
0.9999703305813916
0.9999754797940889
0.9999801386147124
0.9999843070405215
0.9999879850690642
0.9999911726981767
0.9999938699259837
0.9999960767508986
0.9999977931716232
0.9999990191871476
0.9999997547967509
Now we can use the methods given in LocalMethods which eats inner products and returns cov to compute the first column of the within group and accross group covariance matrix.
slice 1 within group covariance to first pixel
11 dot.n. 1, ref, Refn 1, ref1 . covn1dotn2 slice 2 within group covariance to first pixel
22 dot.n. 2, ref, Refn 2, ref1 . covn1dotn2 slice 2 slice 1 across group covariance to first pixel
12 dot.n. 1, ref, Refn 2, ref1 . covn1dotn2 slice 1 slice 2 across group covariance to first pixel
21 dot.n. 2, ref, Refn 1, ref1 . covn1dotn2
10:
10:
8192element ArrayFloat64,1:
748.7111817196152
7
748.708561033119
748.70069908468
748.6875962074554
748.66925295663
748.6456701093434
748.6168486645637
748.5827898430257
748.5434950871988
748.4989660613383
748.4492046516665
748.394212966757
748.3339933380622
748.3339933380622
748.394212966757
748.4492046516665
748.4989660613383
748.5434950871988
748.5827898430257
748.6168486645637
748.6456701093434
748.66925295663
748.6875962074554
748.70069908468
748.708561033119
figurefigsize10,4
subplot1,2,1
plot ref1:200, 111:200, labelslice 1 plot ref1:200, 221:200, , labelslice 2 xlabelazmuth lag
titlewithin slice covariance
legend
subplot1,2,2
plot ref1:200, 211:200, labelslice 1 plot ref1:200, 121:200, , labelslice 1 xlabelazmuth
titleacross slice crosscovariance
legend
11:
8
11: PyObject matplotlib.legend.Legend object at 0x18562293d0 1.6 Spectral analysis
Using the fact that these 4 block matrices are circulant, we take the scaled Discrete Fourier Transform to diagonalize them
12:
12:
Here are the pixel and grid parameters
period 2
period N pixel spacing 2periodFourierspacing nyq
ref . 0:N1 collect
nyq . ifelse ref nyq, ref, ref 2nyq
8192element ArrayFloat64,1:
0.0
1.0
2.0
3.0
4.0
5.0
6.0
7.0
8.0
9.0
10.0
11.0
12.0
9
12.0
11.0
10.0
9.0
8.0
7.0
6.0
5.0
4.0
3.0
2.0
1.0
Here are the fft operators
planfftsimilarcmb1 unscaled fft
1sqrtN
2 sc 2
unitary version
scaled fft which approximates the integral
scaling used for the eigenvalues from the first
col of the circulantmatrix
13:
13:
14:
14:
0.0007669903939428206 FFTW forward plan for 8192element array of
ComplexFloat64
These model cmb1 and cmb2
sc11 sc 11 . real sc22 sc 22 . real sc 21 sc 21
sc 12 sc 12
8192element ArrayComplexFloat64,1:
154.62796786456275 2.786085204295846e16im
599.569757606981 1.6779814513390042e13im
779.6420483508 4.583568012836338e13im
372.8710732243214 3.213635898096056e13im
200.89061308534298 2.432443930829525e13im
116.5589092472053 1.606333851540168e13im
70.9989149940158 1.3088744033495058e13im
44.79815710632035 7.668671487384125e14im
28.98581845967752 6.384661237450354e14im
19.12080169053594 4.618813712709086e14im
12.844497256954483 4.475246054726254e14im
8.730203340362827 3.0643558407250344e14im
5.974482487946962 1.7940551232936005e14im
5.974482487946963 1.9504201137134958e14im
8.730203340362824 3.951261946811799e14im
10
12.844497256954469 3.954397838664373e14im
19.12080169053594 4.8859926922865755e14im
28.98581845967752 6.384412657781292e14im
44.79815710632037 7.424454965801742e14im
70.9989149940158 1.2573664357393976e13im
116.55890924720532 1.6118877847732418e13im
200.890613085343 2.435562912903662e13im
372.8710732243214 3.267591050513847e13im
779.6420483508 4.569026892541495e13im
599.5697576069812 1.681117049520985e13im
The theory of circulant matrices which predicts 2 ij is the expected power in cmbi conj cmbj.
15:
15:
cmb1 cmb1 cmb2 cmb2
8192element ArrayComplexFloat64,1:
133.90814399838476 2.1225283178226174e16im
4.598592707638935 14.63128249815512im
11.703730351827819 19.90828970987386im
37.33998677029809 8.809719862382574im
37.647894786328564 24.40948634772287im
27.78505691144916 13.20704320742245im
3.7853845896430225 6.4362876647001475im
7.8393694627857045 11.794714499160008im
13.234914335385948 15.317386211599176im
1.6906889666819944 12.45987370140233im
8.969810355925363 7.910895730786762im
18.591384350538966 11.662364168115035im
17.818932365135325 23.169844939092187im
17.818932365135325 23.169844939092187im
18.59138435053897 11.662364168115039im
8.969810355925361 7.910895730786761im
1.6906889666819942 12.45987370140233im
13.234914335385948 15.31738621159918im
7.8393694627857045 11.794714499160008im
3.7853845896430203 6.436287664700148im
27.78505691144916 13.207043207422455im
37.647894786328564 24.40948634772287im
37.33998677029809 8.809719862382572im
11.703730351827819 19.90828970987386im
4.598592707638937 14.631282498155121im
We can look at the power in these maps, i.e. cmb1k2 and cmb2k2, compared to what the signal covariance predicts
11
let pltrange 10:N21,
k refpltrange,
km k.2,
abs2 cmb1 abs2. cmb1, abs2 cmb2 abs2. cmb2, sc 11 sc 11,
sc 22 sc 22
formyplot semilogx,loglog
figurefigsize10,4
sub subplot1,2,1
myplotk, km . abs2 cmb1pltrange, labelabs2 cmb1, colorcolors1, ., alpha0.15
colorcolors2
. sc 11pltrange
myplotk, km . sc 11pltrange, labelsc 11, xlabelwavenumber
sub.setylimkm . sc 11pltrange1, 2maximumkm
sub.setxlimk1, maximumk
legend
sub subplot1,2,2
myplotk, km . abs2 cmb2pltrange, labelabs2 cmb2, colorcolors1, ., alpha0.15
myplotk, km . sc 22pltrange, labelsc 22, xlabelwavenumber
sub.setylimkm . sc 22pltrange1, 2maximumkm sub.setxlimk1, maximumk
legend
end for loop over myplot
colorcolors2
. sc 22pltrange
end
16:
12
1.7 Spectral impact of the noise and beam
From https:wiki.cosmos.esa.intplanckpla2015index.phpSummaryofHFIdatacharacteristics
Effective beam solid angle arcmin2 60.44 Effecitive beam FWHM arcmin 7.3 Sensitivity per beam solid angle KCMB 4.3
Temp sensitivity Kdeg
0.55
17:
Compute the standard deviation of white noise in each Healpix pixel. Note:theHealpixpixelarea 1.72arcmin2withNside2048.
13
For the homework … assign the variable sdofnoiseper1p7arcminpixel
the correct value.
Note:theHealpixpixelarea 1.72arcmin2withNside2048.
sdofnoiseper1p7arcminpixel ?????????? bFWHMarcmin 7.3 arcmin
17: 7.3 18:
18:
, let n sdofnoiseper1p7arcminpixel, bfwhm bFWHMarcmin fill n2, size ref
. notethescaling bFWHMrad deg2radbfwhm60
.expbFWHMrad2abs2nyq16log2
,
end
0.2890157703420252, 0.2890157703420252, 0.2890157703420252,
0.2890157703420252, 0.2890157703420252, 0.2890157703420252, 0.2890157703420252,
0.2890157703420252, 0.2890157703420252, 0.2890157703420252 …
0.2890157703420252, 0.2890157703420252, 0.2890157703420252, 0.2890157703420252,
0.2890157703420252, 0.2890157703420252, 0.2890157703420252, 0.2890157703420252,
0.2890157703420252, 0.2890157703420252, 1.0, 0.9999995934139979,
0.9999983736569836, 0.9999963407319328, 0.9999934946438046, 0.9999898353995421,
0.9999853630080723, 0.9999800774803053, 0.9999739788291351, 0.9999670670694394
… 0.999959342218079, 0.9999670670694394, 0.9999739788291351,
0.9999800774803053, 0.9999853630080723, 0.9999898353995421, 0.9999934946438046,
0.9999963407319328, 0.9999983736569836, 0.9999995934139979
Now we can plot the theory beamed signal, the noise and binned bandpowers
let binpowerswidth 12, pltrange 10:N21,
k refpltrange,
km k.2,
abs2 cmb1 abs2. cmb1, abs2 cmb2 abs2. cmb2, sc 11 sc 11,
sc 22 sc 22
19:
14
beamsc 11 sc 11 . abs2. beamsc 22 sc 22 . abs2.
totsc 11 sc 11 . abs2. . totsc 22 sc 22 . abs2. .
formyplot semilogx,loglog
figurefigsize10,4
sub subplot1,2,1
myplotk,
km . abs2 cmb1pltrange xLocalMethods.
binmeanx; bin binpowerswidth,
labelbinned abs2 cmb1, colorcolors1, .,
myplotk, km . pltrange, label N , k, alpha
0.8
myplotk, km . totsc 11pltrange, , labeltotsc 11, colorcolors7
xlabelwavenumber
sub.setylimkm . beamsc 11pltrange1, 2maximumkm . beamsc 11pltrange
sub.setxlimk1, maximumk
legend
sub subplot1,2,2
myplotk,
km . abs2 cmb2pltrange xLocalMethods.
binmeanx; bin binpowerswidth,
labelbinned abs2 cmb2, colorcolors1, .,
myplotk, km . pltrange, label N , k, alpha
0.8
myplotk, km . totsc 22pltrange, , labeltotsc 22, colorcolors7
xlabelwavenumber
sub.setylimkm . beamsc 22pltrange1, 2maximumkm . beamsc 22pltrange
sub.setxlimk1, maximumk
legend
alpha0.15
labelbeamsc 11, colorcolors2
alpha0.15
labelbeamsc 22, colorcolors2
myplotk, km . beamsc 11pltrange,
myplotk, km . beamsc 22pltrange,
15
end
end for loop over myplot
20:
Also the cross correlation. The circulant matrix theory predicts that E cmb1 . conj cmb2 sc12.abs2.andEcmb2.conjcmb1 sc21.abs2.
let binpowerswidth 0, pltrange 2:N21,
k refpltrange,
km k.0,
cross cmb12 cmb1 . conj cmb2, cross cmb21 cmb2 . conj cmb1,
16
sc 12 sc 12, sc 21 sc 21
rcrosscmb crosscmb12.crosscmb21.2.realthismustbe real
icrosscmb crosscmb12.crosscmb21.2.imagthismustbe purely imaginary
rbeamcrosssc real.sc 12 . sc 21 . 2 . abs2. ibeamcrosssc imag.sc 12 . sc 21 . 2 . abs2.
myplot semilogx plot figurefigsize10,4 sub subplot1,2,1
myplotk,
km . rcross cmbpltrange xLocalMethods.
binmeanx; bin binpowerswidth,
labelbinned rcross cmb, colorcolors1, .,
alpha0.15
myplotk, km . rbeamcrosssc pltrange, labelrbeamcrosssc , colorcolors2
xlabelwavenumber
sub.setylimkm . rbeamcrosssc pltrange1, 2maximumkm . rbeamcrosssc pltrange
sub.setxlimk2, maximumk
legend
sub subplot1,2,2
myplotk,
km . icross cmbpltrange xLocalMethods.
binmeanx; bin binpowerswidth,
labelbinned icross cmb, colorcolors1, .,
alpha0.15
myplotk, km . ibeamcrosssc pltrange, labelibeamcrosssc , colorcolors2
xlabelwavenumber
sub.setylimkm . ibeamcrosssc pltrange1, 2maximumkm . ibeamcrosssc pltrange
end
sub.setxlimk2, maximumk
legend
17
20: PyObject matplotlib.legend.Legend object at 0x184e18a490
1.8 Frequencybyfrequency decorrelation of cmb1 and cmb2
21:
totsc 11 sc 11 . abs2. . totsc 22 sc 22 . abs2. . totsc 21 sc 21 . abs2. . real
totsc 21 . sqrt. totsc 11 . totsc 22 nrm cmb1 cmb1 . sqrt.totsc 11
nrm cmb2 cmb2 . sqrt.totsc 22
z1 nrm cmb1 . exp.im.angle. . nrm cmb2 z2 nrm cmb1 . exp.im.angle. . nrm cmb2
1 1 . abs. 2 1 . abs.
let binpowerswidth 4, pltrange 10:N21,
k refpltrange,
abs2z1 abs2.z1,
abs2z2 abs2.z2, crossz2z1 z2 . conj.z1,
1 1, 22
figurefigsize10,4
sub subplot1,2,1
plot
18
binpowerswidth, 15
binpowerswidth, 15
alpha0.15
end
legend
k,
abs2z1pltrange xLocalMethods.binmeanx; bin
labelbinned abs2z1, colorcolors1, ., alpha0. plotk, 1pltrange, label 1, colorcolors2
sub subplot1,2,2
plot
k,
abs2z2pltrange xLocalMethods.binmeanx; bin
labelbinned abs2z2, colorcolors1, ., alpha0. plotk, 2pltrange, label 2, colorcolors2
figurefigsize10,4
plot
k,
real.crossz2z1pltrange xLocalMethods. binmeanx; bin binpowerswidth,
labelrealz1 conjz2, colorcolors1, .,
plot
k,
imag.crossz2z1pltrange xLocalMethods. binmeanx; bin binpowerswidth,
alpha0.15
labelimagz1 conjz2, colorcolors2, .,
plotk, fill0, sizek, k, alpha 0.8
legend
19
21: PyObject matplotlib.legend.Legend object at 0x1856dea910 1.9 Wiener filtered south pole cap
First we multiply the data i.e. the two circular slices stacked by the inverse cov matrix of the data.
Now we make a function that computes the wiener filter on any unobserved circular slice.
22:
For the homework …
Create a function wfcmbfun that takes in a polar angle and
20
returns the conditional expected value of cmb at the azmuth grid values ref.
????????????????
22: wfcmbfun generic function with 1 method Finally we create a grid of polar angles to predict on
23:
if at north pole
s range0.001, max 2, 1 . 0.1, length250 transpose
if at south pole
s rangemin 2, 1 . 0.1, , length150 transpose grid s.0. ref
grid0. s. ref
function makewfcmbmat s, ref wfcmbmatlocal 0 . s . 0 . ref showprogress for icol 1:length s
printicol,
wfcmbmatlocal:,icol wfcmbfun sicol wfcmbmatlocal
end end
wfcmbmat makewfcmbmat s, ref
23:
Progress: 100 Time: 0:06:09
8192150 ArrayFloat64,2:
41.6078 43.7738 45.9311 48.05 50.0846 … 108.516 108.587 108.627
41.9043 44.0884 46.2649 48.4048 50.4616
42.2009 44.4032 46.5991 48.7604 50.8397
42.4976 44.7181 46.9337 49.1168 51.219
42.7944 45.0332 47.2689 49.4739 51.5995
43.0912 45.3486 47.6045 49.8318 51.9812 … 108.515 108.586 108.627
43.3882 45.664 47.9406 50.1903 52.3642
43.6851 45.9797 48.2771 50.5495 52.7483
43.9821 46.2956 48.6141 50.9095 53.1337
44.279 46.6116 48.9516 51.2701 53.5204
44.5759 46.9279 49.2894 51.6315 53.9081 … 108.514 108.586 108.627
108.516 108.587 108.627
108.516 108.587 108.627
108.515 108.587 108.627
108.515 108.587 108.627
108.515 108.586 108.627
108.515 108.586 108.627
108.515 108.586 108.627
108.514 108.586 108.627
21
44.8728 47.2443 49.6277 51.9935 54.297 108.514 108.586 108.627
45.1696 47.5609 49.9663 52.3562 54.6871 108.514 108.586 108.627
38.0591 40.0195 41.9633 43.8583 45.6599 … 108.518 108.588 108.627
38.3541 40.3307 42.2911 44.2028 46.0215
38.6492 40.6422 42.6193 44.5481 46.3844
38.9444 40.9541 42.9482 44.8944 46.7486
39.2398 41.2662 43.2775 45.2415 47.1141
39.5353 41.5787 43.6074 45.5895 47.481 … 108.517 108.588 108.627
39.831 41.8915 43.9378 45.9384 47.8491
40.1268 42.2046 44.2688 46.2882 48.2185
40.4228 42.5179 44.6003 46.6389 48.5892
40.7189 42.8316 44.9323 46.9904 48.9612
41.0151 43.1454 45.2647 47.3428 49.3344 … 108.516 108.587 108.627
41.3114 43.4595 45.5977 47.696 49.7089 108.516 108.587 108.627
108.518 108.588 108.627
108.518 108.588 108.627
108.517 108.588 108.627
108.517 108.588 108.627
108.517 108.587 108.627
108.517 108.587 108.627
108.517 108.587 108.627
108.516 108.587 108.627
fig figure
ax mplot3d.Axes3Dfig
subplotprojectionpolar
if at north pole
pcolormesh grid, grid, wfcmbmat, vmin125, vmax225
if at south pole
pcolormeshgrid, .grid,wfcmbmat,vmin135,vmax235 grid
24:
22
24: PyObject matplotlib.collections.QuadMesh object at 0x14a952b10
23