In [259]:
%matplotlib inline
# We use CAMB generated scalar C_l values, N_side=8
# and CAMB simulated maps associated with scalar C_l values above
# N_side=8 gives N_pix=768
#
# Provide three inputs:
# 1. scalar Cls FITS
# 2. sky map FITS
# 3. manually insert Nside value
# If you did that, alles in Butter!
In [260]:
# http://healpy.readthedocs.org/en/latest/generated/healpy.sphtfunc.anafast.html
#
# lmax = (3*Nside - 1)
#
# if lmax=0, there is only one Cl, i.e. C0. If lmax=11, there are 12 Cl values,
# Cl= [C0, C1 ... C11]
In [261]:
# The standard I am using is lmax = (3*Nside - 1)
#
# However, for all computations, we CONSISTENTLY remove monopole and dipole terms
# Otherwise, we end up dividing by zero
#
# So, generate FITS with lmax = 3Nside+1
In [262]:
import math
import matplotlib.pyplot as plt
import numpy as np
import healpy as hp
import pyfits as pf
import astropy as ap
import os
from scipy.special import eval_legendre ##special scipy function
np.set_printoptions(threshold='nan') # Default is threshold=1000
## Use this to print all values, disables corner printing
In [263]:
cd ~/downloads
In [264]:
print "Changed directory to 'Downloads' folder"
print 'Example of FITS file for C_l scalars is format "camb_cls_nside8_lmax25.fits" '
print 'Example of FITS file for fully sky maps is format "camb_map_nside8_lmax25.fits" '
In [265]:
filename = "camb_cls_nside8_lmax25.fits" # CAMB C_l scalars
temp = "camb_map_nside8_lmax25.fits" # CAMB simulated maps
nside = 8
In [266]:
npix = 12*(nside**2) #total number of pixels, npix
lmax = ((3*nside)-1) #maximum l of the power spectrum C_l
print "The total number of pixels is " + str(npix)
print "The maximum ell of the power spectrum C_l is " + str(lmax)
In [267]:
# open a FITS file
# type()=pyfits.hdu.hdulist.HDUList
ff = pf.open(filename)
In [268]:
### Recall there are four columns: temp, E pol, B pol, grad-temp cross terms
## first two values are zero, i.e. monopole, dipole
cls = ff[1].data # actually second HDU, first is empty
In [269]:
# XXX.field() references columns by 0-index
# field(0) is temperature values
# all Cl scalar temp values put into ndarray
# type()=numpy.ndarray
cltemp = cls.field(0)
newcls = np.delete( cltemp, [0,1] )
# Note: imported C_l values 3*Nside + 1 = cltemp
# newcls uses lmax = 3Nside - 1
print len(cltemp)
print "the number of l values: l=0, l=1, ... l=lmax, where lmax = " + str(lmax)
print "this means length = " + str(len(newcls))
print len(newcls) #the l=0 value
In [270]:
# define ell values
# array from 0 to lmax, where lmax =
# add two values in case there are issues with monopole, dipole term
ellval = np.arange(len(cltemp))
In [271]:
# P_0 is the monopole, P_1 is the dipole
# remove 0, 1
newellval = np.delete(ellval, [0,1]) #numpy.ndarray
print ellval
print newellval
In [272]:
# First calculate the covariance matrix by the definition, i.e.
#
# C_ij = < \delta T_i \delta T_j > = 1/4pi \sum^{N_pix} p =
# ( T_i(p) - mean(T_i) ) * ( T_j(p) - mean(T_j) )
#
#
# Larson, Weiland, Hinshaw, Bennett,
# http://arxiv.org/pdf/1409.7718.pdf
#
In [273]:
tempmap = hp.mrdfits(temp) #type() = list
In [274]:
tempval = tempmap[0][:] #len()=12288, type()=numpy.ndarray
# len(tempval) = 12
# 1024 * 12 = 12288
# 12 rows of 1024 entries
# tempval.shape = (12, 1024)
In [275]:
# The HEALPix standard is fixed at shape = (n, 1024)
# Columns are consistently 1024
# Rows vary
#
# Concatenate array so it is one-dimensional
#
#
# Could use numpy.reshape
# Quicker, probably
# np.reshape(tempval, (1, -1))
#
if npix < 1024:
tempdata = tempval
else:
tempdata = np.concatenate([tempval[i,:] for i in range(tempval.shape[0])])
In [276]:
tempdata.shape
Out[276]:
In [277]:
### Begin calculating S_ij piece by piece, in order to do the summation correctly
#
# S_ij = sum(2ell+1) C_l P_l(dotproductmatrix)
In [278]:
totalpix = np.arange(npix)
#print ell = [ 0 1 2 ..., 3071, 3072]
In [279]:
## healpy.pixelfunc.pix2vec(nside, ipix, nest=False)
##
## will give three arrays
## arrays of all x values, all y values, all z values
## RING scheme default
# len()=3
# type()=tuple
#
#
vecval = hp.pix2vec(nside, totalpix) #Nside = 16, type()=tuple
In [280]:
vecvalx = vecval[0] #shape (3072,)
vecvaly = vecval[1]
vecvalz = vecval[2]
In [281]:
## First arrange arrays vertically
## numpy.vstack = Stack arrays in sequence vertically (row wise), input sequence of arrays
totalvecval = np.vstack((vecvalx, vecvaly, vecvalz)) #type()=numpy.ndarray
In [282]:
trans = totalvecval.T #transpose
In [283]:
dotproductmatrix = trans.dot(totalvecval) #take the dot product
#dotproductmatrix.shape = (npix, npix)
In [284]:
## Begin calculating S_ij piece by piece, in order to do the summation correctly.
#
# S_ij = sum(2ell+1) C_l P_l(dotproductmatrix)
# NOT QUICK!
summatrix = np.sum( eval_legendre(newellval[:, None, None], dotproductmatrix), axis=0)
#
# numpy.sum(a, axis=None, dtype=None, out=None, keepdims=False)[source]
# a : array_like
# Elements to sum.
# axis : None or int or tuple of ints, optional
# Axis or axes along which a sum is performed.
# The default (axis = None) is perform a sum over all the dimensions of the input array.
# axis may be negative, in which case it counts from the last to the first axis.
#
print summatrix.shape
In [285]:
# matrix_total =
# (1/(4*math.pi)) * sum((2 * ll + 1) * cltemp ) * eval_legendre(ll, matrix_dotprod)
#
# Begin with adding theoretical scalar C_l values
#
add_clvalues = np.sum(newcls[:, None, None] * summatrix, axis=0)
print add_clvalues.shape
In [286]:
# matrix_total =
# (1/(4*math.pi)) * np.sum((2*ll + 1) * cltemp ) * eval_legendre(ll, matrix_dotprod)
wholematrix = np.sum( (((2 * newellval[:, None, None]) + 1) * add_clvalues), axis=0)
print wholematrix.shape
In [287]:
covmatrix = (1/(4 * math.pi)) * wholematrix #covariance matrix for Nside=4
In [288]:
# We now compute likelihood functions
#
# For equations, refer to Nside4_covariance_likelihood.tex
# "Nside = 4: Covariance and likelihood"
#
In [ ]:
In [289]:
map1 = hp.synfast(newcls, nside)
map2 = hp.synfast(cltemp, nside)
print map1.shape
tempentries = map2
### healpy.sphtfunc.synfast(cls, nside, lmax=None, mmax=None, alm=False, pol=True,
### pixwin=False, fwhm=0.0, sigma=None, new=False, verbose=True)
# INPUT cls, nside
# OUTPUT sky maps
viewmap = hp.mollview(map1) ##this map DOES NOT contains monopole, dipole
viewmap2 = hp.mollview(map2) ##this map DOES INDEED contain monopol, dipole
In [290]:
powerspectrum = hp.anafast(map1) # create powerspectrum \hat{C}_l
pwer2 = hp.anafast(map2)
print len(powerspectrum)
print powerspectrum
print len(pwer2)
print pwer2
print "***"
print "Notice the first two values, power e-18 and e-16, monopole and dipole"
In [300]:
print "Evaluate spherical harmonic likelihood function"
print "-2\ln L = \sum_{l}(2l+1)\Big[\ln\Big( \frac{C^{\text{th}}_l}{\hat{C}_l}\Big) +\Big(\hat{C}_l/C^{\text{th}}_l\Big)-1 \Big]"
print "where HEALPix anafast is \hat(C)_l, frac{1}{2l+1}\sum_m \vert\hat{a}_{lm}\vert^2"
print "****"
print "We foolishly guess-timate Sum_ell(2*ell +1)"
print "Calculate sum( [(2*i + 1) for i in ell])"
estimatesum = np.sum( [(2*i + 1) for i in newellval], axis=0) #newellval HAS NO monopole, dipole
print "Our crazy guess is something around the value = " + str(estimatesum)
#C^th_l is newcls, NO MONOPOLE, NO DIPOLE
#\hat(C)_l is powerspectrum = hp.anafast(map1), where map1 = hp.synfast(newcls, nside)
#type(newcls) = numpy.ndarray
#type(powerspectrum) = numpy.ndarray
spherlikepart = ( (np.log(newcls/powerspectrum)) + (powerspectrum/newcls) - 1 )
print spherlikepart.shape
sphericalharmoniclikelihood = np.sum( [(2*i +1 ) * spherlikepart for i in newellval])
print sphericalharmoniclikelihood
print "**************"
print "**************"
In [292]:
# Compute temperature space (real space) likelihood
#
# where T is a temperature map array (3072,), S is our above matrix = entirething,
# and N is the number of pixels in a vector
#
# First Year Wilkinson Microwave Anisotropy Probe (WMAP) Observations:
# Parameter Estimation Methodology L. Verde, et al., 2003, ApJS, 148, 195
#
# http://lambda.gsfc.nasa.gov/product/map/dr1/
# pub_papers/firstyear/methodology/wmap_param_method.pdf
#
# print "−2 ln L \propto T*S^{-1}*T + ln detS + N ln 2\pi"
In [293]:
tempvalues = np.matrix(tempentries) #create matrix of temperature values from CAMB map
# tempvalues.shape = (1, 192)
print tempvalues.shape
In [294]:
temptranspose = tempvalues.T
# temptranspose.shape = (192, 1)
In [295]:
## Invert covariance matrix, S_ij
## numpy.linalg.inv(a)
## Compute the (multiplicative) inverse of a matrix.
## Given a square matrix a, return the matrix ainv satisfying
## dot(a, ainv) = dot(ainv, a) = eye(a.shape[0]).
##
inverse = np.linalg.inv(covmatrix)
inversematrix = np.matrix(inverse) # S^{-1}
inversematrix.shape
Out[295]:
In [296]:
# WARNING: DO NOT USE NP.LIGALG.DET(A)!
# Operatioin will not succeed, gives inf or zero
#
#
#
# numpy.linalg.slogdet
# numpy.linalg.slogdet(a)[source]
# Compute the sign and (natural) logarithm of the determinant of an array.
# INPUT array
# OUTPUT (sign, ln det a)
# sign = A number representing the sign of the determinant.
# For a real matrix, this is 1, 0, or -1.
# For a complex matrix, this is a complex number with absolute value 1 (i.e., it is on the unit circle), or else 0.
# ln det a = The natural log of the absolute value of the determinant.
#
lnDetS = np.linalg.slogdet(inversematrix) #type()=tuple
print lnDetS
type(lnDetS[1])
Out[296]:
In [297]:
type(temptranspose)
temptranspose.shape
Out[297]:
In [298]:
realtemplikehd = ( tempvalues * inversematrix * temptranspose )
realtemplikehdpt2 = ( tempvalues * inversematrix * temptranspose ) + lnDetS[1]
realtemplikehdpt3 = ( tempvalues * inversematrix * temptranspose ) + ( lnDetS[1] ) + ((192)*2*math.pi)
In [299]:
print "Real space (temperature space) likelihood function"
print "Formalism: −2 ln L \propto T*S^{-1}*T + ln detS + N ln 2\pi"
print "First number: T*S^(-1)*T "
print realtemplikehd
print "Second number: T*S^{-1}*T + ln detS "
print realtemplikehdpt2
print "Third number: −2 ln L \propto T*S^{-1}*T + ln detS + N ln 2\pi"
print realtemplikehdpt3
print "**************"
print "**************"
print "END"
In [ ]:
In [ ]:
In [ ]:
In [ ]:
In [ ]:
In [ ]: