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.mdfor markdown format convertdocnotebook2.jl, notebook2.ipynbfor 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 datadirUsersethananderesDropboxCoursesSTA250fall2019Datasets Dataset2Planck1dtransectdata
const clsJLD2.FileIO.loaddatadircls.jld2, cls
const wonPlcls:cllenscalar:,1 . 2 . cls:ell . 1 . 4.
cls:factoronclcmb
const lmaxmaximumcls:ell
function covcos cos ,lmax::Int
assert lmax2 lmax too small, must be greater than 1 dotsflegendrePlarraylmax, cos 3:end, wonPl3:lmax1
end
covn1dotn2n1dotn2covcos n1dotn2,lmax
1
covangleanglecovcos cosangle,lmax
function binmeanfk; bin0 fkrtncopyfk
if bin1
Nm1lengthfk
subNm1binNm1bin
fmkreshapefk1:subNm1, bin, Nm1bin fmk . meanfmk, dims1
fkrtn1:subNm1 . vecfmk fkrtnsubNm11:end . meanfksubNm11:end
end
return fkrtn
end
endmodule
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 propcycleplt.rcParamsaxes.propcycle const colorspropcycle.bykeycolor const mplot3dpyimportmpltoolkits.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:
loaditjld2file, varnameJLD2.FileIO.loadLocalMethods. datadirjld2file, varname
northploaditdataslices.jld2,northp
southploaditdataslices.jld2,southp cmbnorthploaditdataslices.jld2, cmbnorthp cmbsouthploaditdataslices.jld2, cmbsouthp
refloaditdataslices.jld2,ref Nlength ref
Nside2048
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:
Islice11
Islice24
1southpIslice1
2southpIslice2 cmb1cmbsouthpIslice1 cmb2cmbsouthpIslice2
show 1, 2
show summarycmb1 show summarycmb2 show summary ref;
1, 21.990924632438072, 2.2330667085254525 summarycmb18192element ArrayFloat64,1 summarycmb28192element ArrayFloat64,1 summary ref8192element ArrayFloat64,1
ofsetcmb9000.vcatnorthp,southp. .0.5reverse figurefigsize10,5
for i1:length northp
plot ref, cmbnorthpi . ofsetcmbi, k, alpha0.1
end
for i1:length southp
if i in Islice1, Islice2
plot ref, cmbsouthpi . ofsetcmbilength southp
else
5:
3
alpha0.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 hppyimporthealpy
const ScSpyimportscipy.spatial const ScIpyimportscipy.interpolate
cmbfullloaditdatafull.jld2, cmbfull fullloaditdatafull.jld2,full fullloaditdatafull.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, ref10.9130366148294711, 0.0, 0.4078776041666667 n 2, ref10.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 wantbroadcasted 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
11dot.n.1, ref, Refn 1, ref1 . covn1dotn2slice 2 within group covariance to first pixel
22dot.n.2, ref, Refn 2, ref1 . covn1dotn2slice 2slice 1 across group covariance to first pixel
12dot.n.1, ref, Refn 2, ref1 . covn1dotn2slice 1slice 2 across group covariance to first pixel
21dot.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
period2
periodNpixel spacing2periodFourierspacing nyq
ref . 0:N1collect
nyq. ifelse refnyq, ref, ref2nyq
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
planfftsimilarcmb1unscaled fft
1sqrtN
2 sc2
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.0007669903939428206FFTW forward plan for 8192element array of
ComplexFloat64
These model cmb1 and cmb2
sc11sc11 . real sc22sc22 . real sc 21sc21
sc 12sc12
8192element ArrayComplexFloat64,1:
154.627967864562752.786085204295846e16im
599.5697576069811.6779814513390042e13im
779.64204835084.583568012836338e13im
372.87107322432143.213635898096056e13im
200.890613085342982.432443930829525e13im
116.55890924720531.606333851540168e13im
70.99891499401581.3088744033495058e13im
44.798157106320357.668671487384125e14im
28.985818459677526.384661237450354e14im
19.120801690535944.618813712709086e14im
12.8444972569544834.475246054726254e14im
8.7302033403628273.0643558407250344e14im
5.9744824879469621.7940551232936005e14im
5.9744824879469631.9504201137134958e14im
8.7302033403628243.951261946811799e14im
10
12.8444972569544693.954397838664373e14im
19.120801690535944.8859926922865755e14im
28.985818459677526.384412657781292e14im
44.798157106320377.424454965801742e14im
70.99891499401581.2573664357393976e13im
116.558909247205321.6118877847732418e13im
200.8906130853432.435562912903662e13im
372.87107322432143.267591050513847e13im
779.64204835084.569026892541495e13im
599.56975760698121.681117049520985e13im
The theory of circulant matrices which predicts 2ij is the expected power in cmbiconj cmbj.
15:
15:
cmb1 cmb1 cmb2 cmb2
8192element ArrayComplexFloat64,1:
133.908143998384762.1225283178226174e16im
4.59859270763893514.63128249815512im
11.70373035182781919.90828970987386im
37.339986770298098.809719862382574im
37.64789478632856424.40948634772287im
27.7850569114491613.20704320742245im
3.78538458964302256.4362876647001475im
7.839369462785704511.794714499160008im
13.23491433538594815.317386211599176im
1.690688966681994412.45987370140233im
8.9698103559253637.910895730786762im
18.59138435053896611.662364168115035im
17.81893236513532523.169844939092187im
17.81893236513532523.169844939092187im
18.5913843505389711.662364168115039im
8.9698103559253617.910895730786761im
1.690688966681994212.45987370140233im
13.23491433538594815.31738621159918im
7.839369462785704511.794714499160008im
3.78538458964302036.436287664700148im
27.7850569114491613.207043207422455im
37.64789478632856424.40948634772287im
37.339986770298098.809719862382572im
11.70373035182781919.90828970987386im
4.59859270763893714.631282498155121im
We can look at the power in these maps, i.e.cmb1k2 andcmb2k2, compared to what the signal covariance predicts
11
let pltrange10:N21,
krefpltrange,
kmk.2,
abs2 cmb1abs2. cmb1, abs2 cmb2abs2. cmb2, sc 11sc 11,
sc 22sc 22
formyplot semilogx,loglog
figurefigsize10,4
subsubplot1,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
subsubplot1,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
endfor 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 anglearcmin2 60.44 Effecitive beam FWHM arcmin 7.3 Sensitivity per beam solid angleKCMB 4.3
Temp sensitivityKdeg
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?????????? bFWHMarcmin7.3arcmin
17: 7.3 18:
18:
,let nsdofnoiseper1p7arcminpixel, bfwhmbFWHMarcminfill n2, size ref
. notethescaling bFWHMraddeg2radbfwhm60
.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 binpowerswidth12, pltrange10:N21,
krefpltrange,
kmk.2,
abs2 cmb1abs2. cmb1, abs2 cmb2abs2. cmb2, sc 11sc 11,
sc 22sc 22
19:
14
beamsc 11sc 11 . abs2.beamsc 22sc 22 . abs2.
totsc 11sc 11 . abs2.. totsc 22sc 22 . abs2..
formyplot semilogx,loglog
figurefigsize10,4
subsubplot1,2,1
myplotk,
km . abs2 cmb1pltrangexLocalMethods.
binmeanx; binbinpowerswidth,
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
subsubplot1,2,2
myplotk,
km . abs2 cmb2pltrangexLocalMethods.
binmeanx; binbinpowerswidth,
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
endfor loop over myplot
20:
Also the cross correlation. The circulant matrix theory predicts that E cmb1 . conj cmb2sc12.abs2.andEcmb2.conjcmb1 sc21.abs2.
let binpowerswidth0, pltrange2:N21,
krefpltrange,
kmk.0,
cross cmb12cmb1 . conj cmb2, cross cmb21cmb2 . conj cmb1,
16
sc 12sc 12, sc 21sc 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.
myplotsemilogx plot figurefigsize10,4 subsubplot1,2,1
myplotk,
km . rcross cmbpltrangexLocalMethods.
binmeanx; binbinpowerswidth,
labelbinned rcross cmb, colorcolors1, .,
alpha0.15
myplotk, km . rbeamcrosssc pltrange, labelrbeamcrosssc , colorcolors2
xlabelwavenumber
sub.setylimkm . rbeamcrosssc pltrange1, 2maximumkm . rbeamcrosssc pltrange
sub.setxlimk2, maximumk
legend
subsubplot1,2,2
myplotk,
km . icross cmbpltrangexLocalMethods.
binmeanx; binbinpowerswidth,
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 11sc 11 . abs2.. totsc 22sc 22 . abs2.. totsc 21sc 21 . abs2.. real
totsc 21 . sqrt. totsc 11 . totsc 22 nrm cmb1cmb1 . sqrt.totsc 11
nrm cmb2cmb2 . sqrt.totsc 22
z1nrm cmb1 . exp.im.angle.. nrm cmb2 z2nrm cmb1 . exp.im.angle.. nrm cmb2
11 . abs.21 . abs.
let binpowerswidth4, pltrange10:N21,
krefpltrange,
abs2z1abs2.z1,
abs2z2abs2.z2, crossz2z1z2 . conj.z1,
11, 22
figurefigsize10,4
subsubplot1,2,1
plot
18
binpowerswidth, 15
binpowerswidth, 15
alpha0.15
end
legend
k,
abs2z1pltrangexLocalMethods.binmeanx; bin
labelbinned abs2z1, colorcolors1, ., alpha0. plotk, 1pltrange, label 1, colorcolors2
subsubplot1,2,2
plot
k,
abs2z2pltrangexLocalMethods.binmeanx; bin
labelbinned abs2z2, colorcolors1, ., alpha0. plotk, 2pltrange, label 2, colorcolors2
figurefigsize10,4
plot
k,
real.crossz2z1pltrangexLocalMethods. binmeanx; binbinpowerswidth,
labelrealz1conjz2, colorcolors1, .,
plot
k,
imag.crossz2z1pltrangexLocalMethods. binmeanx; binbinpowerswidth,
alpha0.15
labelimagz1conjz2, colorcolors2, .,
plotk, fill0, sizek, k, alpha0.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 wfcmbfunthat takes in a polar angle and
20
returns the conditional expected value of cmb at the azmuth grid valuesref.
????????????????
22: wfcmbfun generic function with 1 method Finally we create a grid of polar angles to predict on
23:
if at north pole
srange0.001, max 2, 1 . 0.1, length250transpose
if at south pole
srangemin 2, 1 . 0.1, , length150transpose grid s.0. ref
grid0. s. ref
function makewfcmbmat s, ref wfcmbmatlocal0 . s . 0 . ref showprogress for icol1:length s
printicol,
wfcmbmatlocal:,icolwfcmbfun sicol wfcmbmatlocal
end end
wfcmbmatmakewfcmbmat s, ref
23:
Progress: 100Time: 0:06:09
8192150 ArrayFloat64,2:
41.607843.773845.931148.0550.0846…108.516108.587108.627
41.904344.088446.264948.404850.4616
42.200944.403246.599148.760450.8397
42.497644.718146.933749.116851.219
42.794445.033247.268949.473951.5995
43.091245.348647.604549.831851.9812…108.515108.586108.627
43.388245.664 47.940650.190352.3642
43.685145.979748.277150.549552.7483
43.982146.295648.614150.909553.1337
44.279 46.611648.951651.270153.5204
44.575946.927949.289451.631553.9081…108.514108.586108.627
108.516108.587108.627
108.516108.587108.627
108.515108.587108.627
108.515108.587108.627
108.515108.586108.627
108.515108.586108.627
108.515108.586108.627
108.514108.586108.627
21
44.872847.244349.627751.993554.297108.514108.586108.627
45.169647.560949.966352.356254.6871 108.514108.586108.627
38.059140.019541.963343.858345.6599…108.518108.588108.627
38.354140.330742.291144.202846.0215
38.649240.642242.619344.548146.3844
38.944440.954142.948244.894446.7486
39.239841.266243.277545.241547.1141
39.535341.578743.607445.589547.481 …108.517108.588108.627
39.831 41.891543.937845.938447.8491
40.126842.204644.268846.288248.2185
40.422842.517944.600346.638948.5892
40.718942.831644.932346.990448.9612
41.015143.145445.264747.342849.3344…108.516108.587108.627
41.311443.459545.597747.696 49.7089 108.516108.587108.627
108.518108.588108.627
108.518108.588108.627
108.517108.588108.627
108.517108.588108.627
108.517108.587108.627
108.517108.587108.627
108.517108.587108.627
108.516108.587108.627
figfigure
ax mplot3d.Axes3Dfig
subplotprojectionpolar
if at north pole
pcolormesh grid, grid, wfcmbmat, vmin125, vmax225
if at south pole
pcolormeshgrid, .grid,wfcmbmat,vmin135,vmax235grid
24:
22
24: PyObject matplotlib.collections.QuadMesh object at 0x14a952b10
23
Reviews
There are no reviews yet.