notebook.community

Edit and run



In [15]:

    
from __future__ import (division, print_function, absolute_import)



In [16]:

    
%matplotlib inline
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



In [17]:

    
# Pixel 42 : [ 0.612372, 0.612372, 0.500000 ]
# Pixel 57 : [ 0.783917, 0.523797, 0.333333 ]
# Pixel 58 : [ 0.523797, 0.783917, 0.333333 ]
# Pixel 74 : [ 0.697217, 0.697217, 0.166667 ]



In [18]:

    
#
# We define this pixel patch
#
# Note: We have to put into IDL format for 3D vectors, i.e. 
# HDIL> query_polygon, 512L, [[0.612372, 0.783917, 0.523797, 0.697217], 
#       [0.612372, 0.523797, 0.783917, 0.697217], [0.500000, 0.333333, 0.333333, 0.166667]], listpix3, nlist3
#



In [19]:

    
#
# Now, save IDL .sav file of listpix3
# Import into Python and run
#



In [20]:

    
# http://docs.scipy.org/doc/scipy-0.16.0/reference/generated/scipy.io.readsav.html
# http://www.astrobetter.com/blog/2009/11/24/read-idl-save-files-into-python/



In [21]:

    
import scipy



In [22]:

    
#
# scipy.io.readsav
#
# scipy.io.readsav(file_name, idict=None, python_dict=False, uncompressed_file_name=None, verbose=False)[source]
#
# Read an IDL .sav file
#
#



In [23]:

    
cd ~/Downloads









    



/Users/evanbiederstedt/Downloads



In [24]:

    
import scipy.io



In [25]:

    
patch_file = scipy.io.readsav('listpix_patch3.sav')



In [ ]:



In [26]:

    
type(patch_file)









    Out[26]:





scipy.io.idl.AttrDict



In [27]:

    
arr3 = patch_file['listpix_patch3']
#print(arr3)



In [28]:

    
type(arr3)









    Out[28]:





numpy.ndarray



In [29]:

    
print(len(arr3)) # pixels total 12476



In [30]:

    
camb_map512 = "camb_map_nside512.fits"



In [31]:

    
camb_map512









    Out[31]:





'camb_map_nside512.fits'



In [32]:

    
nside=512
npix = 12*(nside**2) #total number of pixels, npix
LMAX = ((2*nside)) #maximum l of the power spectrum C_l
heal_npix = hp.nside2npix(nside) # Healpix calculated npix

print("The total number of pixels is " + str(npix))
print("The maximum ell of the power spectrum C_l set to lmax = 2*nside " +str(LMAX))
print("Healpix tells me total number of pixels npix is equal to " + str(heal_npix))









    



The total number of pixels is 3145728
The maximum ell of the power spectrum C_l set to lmax = 2*nside 1024
Healpix tells me total number of pixels npix is equal to 3145728



In [33]:

    
mapread_camb512 = hp.read_map(camb_map512)
#hp.mollview(mapread_camb512)









    



NSIDE = 512
ORDERING = RING in fits file



In [34]:

    
# rename array for convenience
tempval = mapread_camb512
#print tempval

# Data:
#     tempval      # the array of pixel values, (3145728,)



In [35]:

    
print(len(tempval))
print(tempval.shape)



In [36]:

    
#
# We only wish to use the pixels defined in our patch
# These pixel indices are listed in arr3 such that total number pixels total 12476
#
# arr3: this defines pixel indices within patch
# 
# To access pixel indices within array of CMB pixels, just use tempval[arr3]
#
patch=tempval[arr3]



In [37]:

    
print(len(patch))



In [38]:

    
# The log-likelihood
#
# -2lnL \propto m^T C^-1 m + ln det C + N ln (2pi)
#
# First term, m^T C^-1 m is the "model fit term"
# Second term, lndetC is the "complexity penalty"
# Third term, N ln 2pi, a constant
#
# m = tempval
# C = Sij



In [39]:

    
m = patch



In [40]:

    
# Next, create the matrix, n_i /cdot n_j
# solely using Healpy routines, i.e. taking the dot product of the vectors
# The result is "dotproductmatrix"



In [41]:

    
npix









    Out[41]:





3145728



In [42]:

    
nside









    Out[42]:





512



In [43]:

    
## 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



In [44]:

    
vecval = hp.pix2vec(nside, arr3) #Nside = 512, type()=tuple



In [45]:

    
len(vecval)









    Out[45]:





3



In [46]:

    
vecvalx = vecval[0] #len() = 12476
vecvaly = vecval[1]
vecvalz = vecval[2]



In [47]:

    
# 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 [48]:

    
trans = totalvecval.T #transpose



In [49]:

    
dotproductmatrix = trans.dot(totalvecval) #take the dot product
print(dotproductmatrix.shape) # = (npix, npix) = (12476, 12476)
# type(dotproductmatrix) = np.ndarray



In [50]:

    
#
# The following procedure is for the angular power spectrum, C^th_ell
# However, we are using some cosmological parameter, /alpha
#
#
# =========================================================
# =========================================================
#
# \Sum_l (2*l + 1)/4pi C^th_l P_l (dotproductmatrix)
# sum from l=2 to l=lmax
#
# arrays l = [2 3 4 .... lmax]
#        C_l = [C_2 C_3 .... C_lmax]
#
# The correct way to do the summation:
# 
# Step 1: calculate the matrix
#            M = dotproductmatrix
#
# Step 2: evaluate the function P_l(x) for each entry of the matrix
#         OUTPUT: [P_2(M) P_3(M) P_4(M) .... P_lmax(M) ]
#
# Step 3: (2*l +1)/4pi from l=2 to l=lmax
#          [5/4pi 7/4pi 9/4pi 11/4pi .... 65/4pi ]
#
# Step 4: multiply 
#         [5/4pi*P_2(M) + 7/4pi*P_3(M) +...... + 65/4pi*P_32(M)]
#
#
# Step 5: multiply by theoretical CAMB values, [C_2 C_3    C_31 C_32]
#         [5/4pi**C_2* P_2(M) + 7/4pi*C_3* P_3(M) +...... + 65/4pi*C_32* P_32(M)]
#
# Step 6: This is an array of S_ij for each theory C_l, l=2 to l=32
#         
#
#
# =========================================================
# =========================================================



In [51]:

    
#print(dotproductmatrix)



In [52]:

    
#
# Let's first just take l_max = nside
# so, that's lmax = 512
#



In [53]:

    
# For lmax = 512, we must create an array of ell values, i.e. [0 1 2 3....31 32]
ell = np.arange(513)
#print(ell)
# 
# Subtract the monopole and dipole, l=0, l=1
ellval = ell[2:]
#print(ellval)



In [54]:

    
# The correct way to do the summation:
# 
# Step 1: calculate the matrix
#            M = dotproductmatrix
#
# Step 2: evaluate the function P_l(x) for each entry of the matrix
#         OUTPUT: [P_2(M) P_3(M) P_4(M) .... P_lmax(M) ]
#
# Step 3: (2*l +1)/4pi from l=2 to l=lmax
#          [5/4pi 7/4pi 9/4pi 11/4pi .... 65/4pi ]
#
# Step 4: multiply 
#         [5/4pi*P_2(M) + 7/4pi*P_3(M) +...... + 65/4pi*P_32(M)]
#
#
# Step 5: multiply by theoretical CAMB values, [C_2 C_3    C_31 C_32]
#         [5/4pi**C_2* P_2(M) + 7/4pi*C_3* P_3(M) +...... + 65/4pi*C_32* P_32(M)]
#
# Step 6: This is an array of S_ij for each theory C_l, l=2 to l=32
#



In [55]:

    
dotproductmatrix.shape









    Out[55]:





(768, 768)



In [56]:

    
# Step 1: calculate the matrix
M = dotproductmatrix



In [57]:

    
# Step 2: evaluate the function P_l(x) for each entry of the matrix
#         OUTPUT: [P_2(M) P_3(M) P_4(M) .... P_lmax(M) ]



In [58]:

    
# CODE BOTTLENECK!
# 
# Evaluate Legendre from l=2 to l=lmax for each matrix entry
# [P_2(M) P_3(M) P_4(M) .... P_lmax(M) ]
# 
# WITHOUT BROADCASTING, one would do something like 
# PlMat = []
# for i in ellval:
#    PlMat.append( eval_legendre(i, dotproductmatrix) )
#
#
# With broadcasting, we use
# PlMat = eval_legendre(ellval[:, None, None], dotproductmatrix)
# PlMat = [P_2(M) P_3(M) P_4(M) .... P_lmax(M) ]
# PlMat is an array, len()=31 of 31 3072 by 3072 matrices
# PlMat.shape = (31, 3072, 3072)



In [59]:

    
#This doesn't run for lmax=512
#So, split 'ellval' into ten arrays and then sum afterwards

splitell = np.array_split(ellval, 150)
splitell[0]









    Out[59]:





array([2, 3, 4, 5])



In [60]:

    
PlMat1 = eval_legendre(splitell[0][:, None, None], dotproductmatrix)



In [61]:

    
PlMat2 = eval_legendre(splitell[1][:, None, None], dotproductmatrix)



In [62]:

    
PlMat3 = eval_legendre(splitell[2][:, None, None], dotproductmatrix)



In [63]:

    
PlMat4 = eval_legendre(splitell[3][:, None, None], dotproductmatrix)



In [64]:

    
PlMat5 = eval_legendre(splitell[4][:, None, None], dotproductmatrix)



In [65]:

    
PlMat6 = eval_legendre(splitell[5][:, None, None], dotproductmatrix)



In [66]:

    
PlMat7 = eval_legendre(splitell[6][:, None, None], dotproductmatrix)



In [67]:

    
PlMat8 = eval_legendre(splitell[7][:, None, None], dotproductmatrix)



In [68]:

    
PlMat9 = eval_legendre(splitell[8][:, None, None], dotproductmatrix)



In [ ]:



In [69]:

    
PlMat10 = eval_legendre(splitell[9][:, None, None], dotproductmatrix)



In [70]:

    
PlMat11 = eval_legendre(splitell[10][:, None, None], dotproductmatrix)



In [71]:

    
PlMat12 = eval_legendre(splitell[11][:, None, None], dotproductmatrix)



In [72]:

    
PlMat13 = eval_legendre(splitell[12][:, None, None], dotproductmatrix)



In [73]:

    
PlMat14 = eval_legendre(splitell[13][:, None, None], dotproductmatrix)



In [74]:

    
PlMat15 = eval_legendre(splitell[14][:, None, None], dotproductmatrix)



In [75]:

    
PlMat16 = eval_legendre(splitell[15][:, None, None], dotproductmatrix)



In [76]:

    
PlMat17 = eval_legendre(splitell[16][:, None, None], dotproductmatrix)



In [77]:

    
PlMat18 = eval_legendre(splitell[17][:, None, None], dotproductmatrix)



In [78]:

    
PlMat19 = eval_legendre(splitell[18][:, None, None], dotproductmatrix)



In [79]:

    
PlMat20 = eval_legendre(splitell[19][:, None, None], dotproductmatrix)



In [80]:

    
PlMat21 = eval_legendre(splitell[20][:, None, None], dotproductmatrix)



In [81]:

    
PlMat22 = eval_legendre(splitell[21][:, None, None], dotproductmatrix)



In [82]:

    
PlMat23 = eval_legendre(splitell[22][:, None, None], dotproductmatrix)



In [83]:

    
PlMat24 = eval_legendre(splitell[23][:, None, None], dotproductmatrix)



In [84]:

    
PlMat25 = eval_legendre(splitell[24][:, None, None], dotproductmatrix)



In [85]:

    
PlMat26 = eval_legendre(splitell[25][:, None, None], dotproductmatrix)



In [86]:

    
PlMat27 = eval_legendre(splitell[26][:, None, None], dotproductmatrix)



In [87]:

    
PlMat28 = eval_legendre(splitell[27][:, None, None], dotproductmatrix)



In [88]:

    
PlMat29 = eval_legendre(splitell[28][:, None, None], dotproductmatrix)



In [89]:

    
PlMat30 = eval_legendre(splitell[29][:, None, None], dotproductmatrix)



In [90]:

    
PlMat31 = eval_legendre(splitell[30][:, None, None], dotproductmatrix)



In [91]:

    
PlMat32 = eval_legendre(splitell[31][:, None, None], dotproductmatrix)



In [92]:

    
PlMat33 = eval_legendre(splitell[32][:, None, None], dotproductmatrix)



In [93]:

    
PlMat34 = eval_legendre(splitell[33][:, None, None], dotproductmatrix)



In [94]:

    
PlMat35 = eval_legendre(splitell[34][:, None, None], dotproductmatrix)



In [95]:

    
PlMat36 = eval_legendre(splitell[35][:, None, None], dotproductmatrix)



In [96]:

    
PlMat37 = eval_legendre(splitell[36][:, None, None], dotproductmatrix)



In [97]:

    
PlMat38 = eval_legendre(splitell[37][:, None, None], dotproductmatrix)



In [98]:

    
PlMat39 = eval_legendre(splitell[38][:, None, None], dotproductmatrix)



In [99]:

    
PlMat40 = eval_legendre(splitell[39][:, None, None], dotproductmatrix)



In [100]:

    
PlMat41 = eval_legendre(splitell[40][:, None, None], dotproductmatrix)



In [101]:

    
PlMat42 = eval_legendre(splitell[41][:, None, None], dotproductmatrix)



In [102]:

    
PlMat43 = eval_legendre(splitell[42][:, None, None], dotproductmatrix)



In [103]:

    
PlMat44 = eval_legendre(splitell[43][:, None, None], dotproductmatrix)



In [104]:

    
PlMat45 = eval_legendre(splitell[44][:, None, None], dotproductmatrix)



In [105]:

    
PlMat46 = eval_legendre(splitell[45][:, None, None], dotproductmatrix)



In [106]:

    
PlMat47 = eval_legendre(splitell[46][:, None, None], dotproductmatrix)



In [107]:

    
PlMat48 = eval_legendre(splitell[47][:, None, None], dotproductmatrix)



In [108]:

    
PlMat49 = eval_legendre(splitell[48][:, None, None], dotproductmatrix)



In [109]:

    
PlMat50 = eval_legendre(splitell[49][:, None, None], dotproductmatrix)



In [110]:

    
PlMat51 = eval_legendre(splitell[50][:, None, None], dotproductmatrix)



In [111]:

    
PlMat52 = eval_legendre(splitell[51][:, None, None], dotproductmatrix)



In [112]:

    
PlMat53 = eval_legendre(splitell[52][:, None, None], dotproductmatrix)



In [113]:

    
PlMat54 = eval_legendre(splitell[53][:, None, None], dotproductmatrix)



In [114]:

    
PlMat55 = eval_legendre(splitell[54][:, None, None], dotproductmatrix)



In [115]:

    
PlMat56 = eval_legendre(splitell[55][:, None, None], dotproductmatrix)



In [116]:

    
PlMat57 = eval_legendre(splitell[56][:, None, None], dotproductmatrix)



In [117]:

    
PlMat58 = eval_legendre(splitell[57][:, None, None], dotproductmatrix)



In [118]:

    
PlMat59 = eval_legendre(splitell[58][:, None, None], dotproductmatrix)



In [119]:

    
PlMat60 = eval_legendre(splitell[59][:, None, None], dotproductmatrix)



In [120]:

    
PlMat61 = eval_legendre(splitell[60][:, None, None], dotproductmatrix)



In [121]:

    
PlMat62 = eval_legendre(splitell[61][:, None, None], dotproductmatrix)



In [122]:

    
PlMat63 = eval_legendre(splitell[62][:, None, None], dotproductmatrix)



In [123]:

    
PlMat64 = eval_legendre(splitell[63][:, None, None], dotproductmatrix)



In [124]:

    
PlMat65 = eval_legendre(splitell[64][:, None, None], dotproductmatrix)



In [125]:

    
PlMat66 = eval_legendre(splitell[65][:, None, None], dotproductmatrix)



In [126]:

    
PlMat67 = eval_legendre(splitell[66][:, None, None], dotproductmatrix)



In [127]:

    
PlMat68 = eval_legendre(splitell[67][:, None, None], dotproductmatrix)



In [128]:

    
PlMat69 = eval_legendre(splitell[68][:, None, None], dotproductmatrix)



In [129]:

    
PlMat70 = eval_legendre(splitell[69][:, None, None], dotproductmatrix)



In [130]:

    
PlMat71 = eval_legendre(splitell[70][:, None, None], dotproductmatrix)



In [131]:

    
PlMat72 = eval_legendre(splitell[71][:, None, None], dotproductmatrix)



In [132]:

    
PlMat73 = eval_legendre(splitell[72][:, None, None], dotproductmatrix)



In [133]:

    
PlMat74 = eval_legendre(splitell[73][:, None, None], dotproductmatrix)



In [134]:

    
PlMat75 = eval_legendre(splitell[74][:, None, None], dotproductmatrix)



In [135]:

    
PlMat76 = eval_legendre(splitell[75][:, None, None], dotproductmatrix)



In [136]:

    
PlMat77 = eval_legendre(splitell[76][:, None, None], dotproductmatrix)



In [137]:

    
PlMat78 = eval_legendre(splitell[77][:, None, None], dotproductmatrix)



In [138]:

    
PlMat79 = eval_legendre(splitell[78][:, None, None], dotproductmatrix)



In [139]:

    
PlMat80 = eval_legendre(splitell[79][:, None, None], dotproductmatrix)



In [140]:

    
PlMat81 = eval_legendre(splitell[80][:, None, None], dotproductmatrix)



In [141]:

    
PlMat82 = eval_legendre(splitell[81][:, None, None], dotproductmatrix)



In [142]:

    
PlMat83 = eval_legendre(splitell[82][:, None, None], dotproductmatrix)



In [143]:

    
PlMat84 = eval_legendre(splitell[83][:, None, None], dotproductmatrix)



In [144]:

    
PlMat85 = eval_legendre(splitell[84][:, None, None], dotproductmatrix)



In [145]:

    
PlMat86 = eval_legendre(splitell[85][:, None, None], dotproductmatrix)



In [146]:

    
PlMat87 = eval_legendre(splitell[86][:, None, None], dotproductmatrix)



In [147]:

    
PlMat88 = eval_legendre(splitell[87][:, None, None], dotproductmatrix)



In [148]:

    
PlMat89 = eval_legendre(splitell[88][:, None, None], dotproductmatrix)



In [149]:

    
PlMat90 = eval_legendre(splitell[89][:, None, None], dotproductmatrix)



In [150]:

    
PlMat91 = eval_legendre(splitell[90][:, None, None], dotproductmatrix)



In [151]:

    
PlMat92 = eval_legendre(splitell[91][:, None, None], dotproductmatrix)



In [152]:

    
PlMat93 = eval_legendre(splitell[92][:, None, None], dotproductmatrix)



In [153]:

    
PlMat94 = eval_legendre(splitell[93][:, None, None], dotproductmatrix)



In [154]:

    
PlMat95 = eval_legendre(splitell[94][:, None, None], dotproductmatrix)



In [155]:

    
PlMat96 = eval_legendre(splitell[95][:, None, None], dotproductmatrix)



In [156]:

    
PlMat97 = eval_legendre(splitell[96][:, None, None], dotproductmatrix)



In [157]:

    
PlMat98 = eval_legendre(splitell[97][:, None, None], dotproductmatrix)



In [158]:

    
PlMat99 = eval_legendre(splitell[98][:, None, None], dotproductmatrix)



In [159]:

    
PlMat100 = eval_legendre(splitell[99][:, None, None], dotproductmatrix)



In [160]:

    
PlMat101 = eval_legendre(splitell[100][:, None, None], dotproductmatrix)



In [161]:

    
PlMat102 = eval_legendre(splitell[101][:, None, None], dotproductmatrix)



In [162]:

    
PlMat103 = eval_legendre(splitell[102][:, None, None], dotproductmatrix)



In [163]:

    
PlMat104 = eval_legendre(splitell[103][:, None, None], dotproductmatrix)



In [164]:

    
PlMat105 = eval_legendre(splitell[104][:, None, None], dotproductmatrix)



In [165]:

    
PlMat106 = eval_legendre(splitell[105][:, None, None], dotproductmatrix)



In [166]:

    
PlMat107 = eval_legendre(splitell[106][:, None, None], dotproductmatrix)



In [167]:

    
PlMat108 = eval_legendre(splitell[107][:, None, None], dotproductmatrix)



In [168]:

    
PlMat109 = eval_legendre(splitell[108][:, None, None], dotproductmatrix)



In [169]:

    
PlMat110 = eval_legendre(splitell[109][:, None, None], dotproductmatrix)



In [170]:

    
PlMat111 = eval_legendre(splitell[110][:, None, None], dotproductmatrix)



In [171]:

    
PlMat112 = eval_legendre(splitell[111][:, None, None], dotproductmatrix)



In [172]:

    
PlMat113 = eval_legendre(splitell[112][:, None, None], dotproductmatrix)



In [173]:

    
PlMat114 = eval_legendre(splitell[113][:, None, None], dotproductmatrix)



In [174]:

    
PlMat115 = eval_legendre(splitell[114][:, None, None], dotproductmatrix)



In [175]:

    
PlMat116 = eval_legendre(splitell[115][:, None, None], dotproductmatrix)



In [176]:

    
PlMat117 = eval_legendre(splitell[116][:, None, None], dotproductmatrix)



In [177]:

    
PlMat118 = eval_legendre(splitell[117][:, None, None], dotproductmatrix)



In [178]:

    
PlMat119 = eval_legendre(splitell[118][:, None, None], dotproductmatrix)



In [179]:

    
PlMat120 = eval_legendre(splitell[119][:, None, None], dotproductmatrix)



In [180]:

    
PlMat121 = eval_legendre(splitell[120][:, None, None], dotproductmatrix)



In [181]:

    
PlMat122 = eval_legendre(splitell[121][:, None, None], dotproductmatrix)



In [182]:

    
PlMat123 = eval_legendre(splitell[122][:, None, None], dotproductmatrix)



In [183]:

    
PlMat124 = eval_legendre(splitell[123][:, None, None], dotproductmatrix)



In [184]:

    
PlMat125 = eval_legendre(splitell[124][:, None, None], dotproductmatrix)



In [185]:

    
PlMat126 = eval_legendre(splitell[125][:, None, None], dotproductmatrix)



In [186]:

    
PlMat127 = eval_legendre(splitell[126][:, None, None], dotproductmatrix)



In [187]:

    
PlMat128 = eval_legendre(splitell[127][:, None, None], dotproductmatrix)



In [188]:

    
PlMat129 = eval_legendre(splitell[128][:, None, None], dotproductmatrix)



In [189]:

    
PlMat130 = eval_legendre(splitell[129][:, None, None], dotproductmatrix)



In [190]:

    
PlMat131 = eval_legendre(splitell[130][:, None, None], dotproductmatrix)



In [191]:

    
PlMat132 = eval_legendre(splitell[131][:, None, None], dotproductmatrix)



In [192]:

    
PlMat133 = eval_legendre(splitell[132][:, None, None], dotproductmatrix)



In [193]:

    
PlMat134 = eval_legendre(splitell[133][:, None, None], dotproductmatrix)



In [194]:

    
PlMat135 = eval_legendre(splitell[134][:, None, None], dotproductmatrix)



In [195]:

    
PlMat136 = eval_legendre(splitell[135][:, None, None], dotproductmatrix)



In [196]:

    
PlMat137 = eval_legendre(splitell[136][:, None, None], dotproductmatrix)



In [197]:

    
PlMat138 = eval_legendre(splitell[137][:, None, None], dotproductmatrix)



In [198]:

    
PlMat139 = eval_legendre(splitell[138][:, None, None], dotproductmatrix)



In [199]:

    
PlMat140 = eval_legendre(splitell[139][:, None, None], dotproductmatrix)



In [200]:

    
PlMat141 = eval_legendre(splitell[140][:, None, None], dotproductmatrix)



In [201]:

    
PlMat142 = eval_legendre(splitell[141][:, None, None], dotproductmatrix)



In [202]:

    
PlMat143 = eval_legendre(splitell[142][:, None, None], dotproductmatrix)



In [203]:

    
PlMat144 = eval_legendre(splitell[143][:, None, None], dotproductmatrix)



In [204]:

    
PlMat145 = eval_legendre(splitell[144][:, None, None], dotproductmatrix)



In [205]:

    
PlMat146 = eval_legendre(splitell[145][:, None, None], dotproductmatrix)



In [206]:

    
PlMat147 = eval_legendre(splitell[146][:, None, None], dotproductmatrix)



In [207]:

    
PlMat148 = eval_legendre(splitell[147][:, None, None], dotproductmatrix)



In [208]:

    
PlMat149 = eval_legendre(splitell[148][:, None, None], dotproductmatrix)



In [209]:

    
PlMat150 = eval_legendre(splitell[149][:, None, None], dotproductmatrix)



In [210]:

    
splitell[49]









    Out[210]:





array([198, 199, 200, 201])



In [211]:

    
PlMat_total = np.concatenate((PlMat1, PlMat2, PlMat3, PlMat4, PlMat5, PlMat6, PlMat7, 
                                 PlMat8, PlMat9, PlMat10, PlMat11, PlMat12, PlMat13, PlMat14, PlMat15, 
                                PlMat16, PlMat17, PlMat18, PlMat19, PlMat20, PlMat21, PlMat22, PlMat23, 
                                PlMat24, PlMat25, PlMat26, PlMat27, PlMat28, PlMat29, PlMat30, PlMat31, 
                                PlMat32, PlMat33, PlMat34, PlMat35, PlMat36, PlMat37, PlMat38, PlMat39, 
                                 PlMat40, PlMat41, PlMat42, PlMat43, PlMat44, PlMat45, PlMat46, PlMat47,
                                 PlMat48, PlMat49, PlMat50, PlMat51, PlMat52, PlMat53, PlMat54, PlMat55,
                                 PlMat56, PlMat57, PlMat58, PlMat59, PlMat60, PlMat61, PlMat62, PlMat63,
                                 PlMat64, PlMat65, PlMat66, PlMat67, PlMat68, PlMat69, PlMat70, PlMat71,
                                 PlMat72, PlMat73, PlMat74, PlMat75, PlMat76, PlMat77, PlMat78, PlMat79,
                                 PlMat80, PlMat81, PlMat82, PlMat83, PlMat84, PlMat85, PlMat86, PlMat87,
                                 PlMat88, PlMat89, PlMat90, PlMat91, PlMat92, PlMat93, PlMat94, PlMat95,
                                PlMat96, PlMat97, PlMat98, PlMat99, PlMat100, PlMat101, PlMat102, PlMat103,
                                 PlMat104, PlMat105, PlMat106, PlMat107, PlMat108, PlMat109, PlMat110, PlMat111,
                                 PlMat112, PlMat113, PlMat114, PlMat115, PlMat116, PlMat117, PlMat118, PlMat119,
                                 PlMat120, PlMat121, PlMat122, PlMat123, PlMat124, PlMat125, PlMat126, PlMat127,
                                 PlMat128, PlMat129, PlMat130, PlMat131, PlMat132, PlMat133, PlMat134, PlMat135,
                                  PlMat136, PlMat137, PlMat138, PlMat139, PlMat140, PlMat141, PlMat142, PlMat143,
                                 PlMat144, PlMat145, PlMat146, PlMat147, PlMat148, PlMat149, PlMat150))



In [212]:

    
PlMat_total.shape









    Out[212]:





(511, 768, 768)



In [ ]:



In [213]:

    
PlMat = PlMat_total



In [214]:

    
#import cPickle as pickle

#file_Name = "testfileNov18"
# open the file for writing
#fileObject = open(file_Name,'wb') 

# this writes the object a to the
# file named 'testfile'
#pickle.dump(PlMat_total, fileObject) 

# here we close the fileObject
#fileObject.close()



In [215]:

    
# Step 3: (2*l +1)/4pi from l=2 to l=lmax
#          [5/4pi 7/4pi 9/4pi 11/4pi .... 65/4pi ]
norm = ((2*ellval + 1))/(4*math.pi)
print(len(ellval))



In [216]:

    
# Step 4: multiply 
#         [5/4pi*P_2(M) + 7/4pi*P_3(M) +...... + 65/4pi*P_32(M)]
#
# multiply PlMat by (2*l+1)/4pi, i.e. norm
norm_matrix = norm[:, None, None] * PlMat
# [5/4pi * P_2(M)  7/4pi * P_3(M) ....   65/4pi * P_32(M)]



In [217]:

    
norm_matrix.shape









    Out[217]:





(511, 768, 768)



In [218]:

    
PlMat.shape









    Out[218]:





(511, 768, 768)



In [ ]:



In [219]:

    
# Step 5: multiply by theoretical CAMB values, [C_2 C_3    C_31 C_32]
#         [5/4pi**C_2* P_2(M) + 7/4pi*C_3* P_3(M) +...... + 65/4pi*C_32* P_32(M)]

#
# In our case, we use CAMB cosmological parameters
# 
# H_0,  Hubble Constant, 70
#
# Tcmb, 2.7255
#
# Ωbh2, 0.0226
#
# Ωch2, 0.112
#
# Ωνh2, 0.00064, neutrino density
#
# Ωk, 0
#
#The six parameters of the basic ΛCDM model are: the physical baryon density, Ωbh
#2
#; the physical cold dark matter
#density, Ωch
#2
#; the dark energy density, in units of the critical density, ΩΛ; the amplitude of primordial scalar curvature
#perturbations, ∆2
#R at k = 0.002 Mpc−1
#; the power-law spectral index of primordial density (scalar) perturbations,
#ns; and the reionization optical depth, τ . In this model, the Hubble constant, H0 = 100h km/s/Mpc, is implicitly
#determined by the flatness constraint, Ωb + Ωc + ΩΛ = 1. A handful of parameters in this model take assumed values
#that we further test in §4; other parameters may be derived from the fit, as in Table 2. Throughout this paper we
#assume the initial fluctuations are adiabatic and Gaussian distributed (see Bennett et al. (2012) for limits on nonGaussian
#fluctuations from the nine



In [220]:

    
ones_arr = np.ones(511)



In [221]:

    
len(ones_arr)









    Out[221]:





511



In [222]:

    
Hubble_constant = 70*ones_arr
Baryon_density = (0.0226)*ones_arr
CDM_density = (0.112)*ones_arr



In [223]:

    
#
# TEST LogLF vs Cosmological parameters
#
# The covariance matrix is a function of variable "Param", an unknown parameter. 
# 
# Our covariance matrix is therefore S_ij  = (2*l +1)/4pi * Param * P_ell(matrix)   
# 
# 
# The LF is then a function of param, LF(param). This is the only parameter we vary. 
#
# LF = -2loglikelihood /propto T^T inv(S_ij) T + log det (Sij) + N log (2pi)
# 
# We then plot LF(x) vs. parameters x. 
#



In [224]:

    
# define pixel-value arrays
mT = np.matrix(patch)     # mT.shape = (1, 3072)
m = np.matrix(patch).T    # m.shape = (3072, 1)
Npix2pi = (len(patch))*2*math.pi  # LF constant



In [225]:

    
# baryon density is 0.0226



In [226]:

    
vary_baryon = np.linspace(0.0050, 0.050, num=20 ) #set default num = 20



In [227]:

    
vary_baryon









    Out[227]:





array([ 0.005     ,  0.00736842,  0.00973684,  0.01210526,  0.01447368,
        0.01684211,  0.01921053,  0.02157895,  0.02394737,  0.02631579,
        0.02868421,  0.03105263,  0.03342105,  0.03578947,  0.03815789,
        0.04052632,  0.04289474,  0.04526316,  0.04763158,  0.05      ])



In [228]:

    
vary_baryon









    Out[228]:





array([ 0.005     ,  0.00736842,  0.00973684,  0.01210526,  0.01447368,
        0.01684211,  0.01921053,  0.02157895,  0.02394737,  0.02631579,
        0.02868421,  0.03105263,  0.03342105,  0.03578947,  0.03815789,
        0.04052632,  0.04289474,  0.04526316,  0.04763158,  0.05      ])



In [229]:

    
sigma2 = np.logspace(-12, -14, num=20) #set default num = 30



In [230]:

    
print(sigma2)









    



[  1.00000000e-12   7.84759970e-13   6.15848211e-13   4.83293024e-13
   3.79269019e-13   2.97635144e-13   2.33572147e-13   1.83298071e-13
   1.43844989e-13   1.12883789e-13   8.85866790e-14   6.95192796e-14
   5.45559478e-14   4.28133240e-14   3.35981829e-14   2.63665090e-14
   2.06913808e-14   1.62377674e-14   1.27427499e-14   1.00000000e-14]



In [231]:

    
# For N matrix, set the identity
id_mat = np.identity(768)



In [232]:

    
noiseresult = sigma2[:, None, None] * id_mat[None, :, :]



In [233]:

    
noiseresult.shape









    Out[233]:





(20, 768, 768)



In [234]:

    
correctmatrix = np.sum(norm_matrix, axis=0)
print(norm_matrix.shape)
print(correctmatrix.shape)









    



(511, 768, 768)
(768, 768)



In [235]:

    
tempp = (1e6)*patch # multiply CMB maps by 1e6

def LogLikehood_wNoise_1e12(param, sig):
    # param is our parameter, C_3
    Sij = param[:, None, None] * correctmatrix[None, :, :]
    newSij = (1e12)*Sij   # multiply S_ij by 1e12
    Nij = sig[:, None, None] * id_mat[None, :, :]
    newNij = (1e12)*Nij
    # Format 7/4pi * param * P_3(M) where param is the parameter we vary, C_l
    # Sij.shape = (20, 3072, 3072)
    Cij = newSij + newNij
    #invCij = np.linalg.inv(Cij)
    logdetC = np.linalg.slogdet(Cij)  # returns sign and determinant; use logdetC[1]
    # model_fit_terms = m^T C^-1 m
    #
    # model_fit_terms = np.array([np.dot(tempval.T , np.dot(invCij[i] , tempval) ) 
    # for i in range(invCij.shape[0])])
    #
    model_fit_terms = np.array([np.dot(tempp.T , np.linalg.solve(Cij[i], tempp) ) for i in range(Cij.shape[0]) ]) 
    return model_fit_terms + logdetC[1] + Npix2pi



In [236]:

    
tempp = (1e6)*patch # multiply CMB maps by 1e6

def LogLikehood_wNoise_1e12_Contour(param, sig):
    # param is our parameter, C_3
    Sij = param * correctmatrix[None, :, :]
    newSij = (1e12)*Sij   # multiply S_ij by 1e12
    Nij = sig * id_mat[None, :, :]
    newNij = (1e12)*Nij
    # Format 7/4pi * param * P_3(M) where param is the parameter we vary, C_l
    # Sij.shape = (20, 3072, 3072)
    Cij = newSij + newNij
    #invCij = np.linalg.inv(Cij)
    logdetC = np.linalg.slogdet(Cij)  # returns sign and determinant; use logdetC[1]
    # model_fit_terms = m^T C^-1 m
    #
    # model_fit_terms = np.array([np.dot(tempval.T , np.dot(invCij[i] , tempval) ) 
    # for i in range(invCij.shape[0])])
    #
    model_fit_terms = np.array([np.dot(tempp.T , np.linalg.solve(Cij[i], tempp) ) for i in range(Cij.shape[0]) ]) 
    return model_fit_terms + logdetC[1] + Npix2pi



In [237]:

    
import pylab as pb    
import matplotlib.pyplot as plt 

vary_baryon = np.linspace(0.0050, 0.050, num=20 ) #set default num = 20
sigma2 = np.logspace(-12, -14, num=20) #set default num = 30

xxx = vary_baryon
yyy = sigma2

zzz = np.array([[LogLikehood_wNoise_1e12_Contour(np.asarray(i), np.asarray(j)) for i in xxx] for j in yyy])

zzzreshaped = zzz.reshape(20,20)

plt.figure()
CS = plt.contour(xxx, yyy, zzzreshaped)
plt.clabel(CS, inline=1, fontsize=10)
plt.xlabel("vary_baryon = np.linspace(0.0050, 0.050, num=20 ) ")
plt.ylabel("sigma2 = np.logspace(-12, -14, num=20)")
pb.show()



In [238]:

    
patch[:15]









    Out[238]:





array([ -1.11058935e-05,  -1.10077535e-05,  -1.14851618e-05,
        -1.09077973e-05,  -1.14025888e-05,  -1.18686266e-05,
        -1.08059849e-05,  -1.13181786e-05,  -1.18022408e-05,
        -1.22554611e-05,  -1.07022615e-05,  -1.12318739e-05,
        -1.17340069e-05,  -1.22058727e-05,  -1.26448595e-05])



In [239]:

    
vary_baryon = np.linspace(0.0050, 0.050, num=20 ) #set default num = 20
sigma2 = np.logspace(-5, -7, num=20) #set default num = 30

xxx = vary_baryon
yyy = sigma2

zzz = np.array([[LogLikehood_wNoise_1e12_Contour(np.asarray(i), np.asarray(j)) for i in xxx] for j in yyy])

zzzreshaped = zzz.reshape(20,20)

plt.figure()
CS = plt.contour(xxx, yyy, zzzreshaped)
plt.clabel(CS, inline=1, fontsize=10)
plt.xlabel("vary_baryon = np.linspace(0.0050, 0.050, num=20 ) ")
plt.ylabel("sigma2 = np.logspace(-5, -7, num=20)")
pb.show()



In [240]:

    
vary_baryon = np.linspace(0.0050, 0.050, num=20 ) #set default num = 20
sigma2 = np.logspace(-6, -8, num=20) #set default num = 30

xxx = vary_baryon
yyy = sigma2

zzz = np.array([[LogLikehood_wNoise_1e12_Contour(np.asarray(i), np.asarray(j)) for i in xxx] for j in yyy])

zzzreshaped = zzz.reshape(20,20)

plt.figure()
CS = plt.contour(xxx, yyy, zzzreshaped)
plt.clabel(CS, inline=1, fontsize=10)
plt.xlabel("vary_baryon = np.linspace(0.0050, 0.050, num=20 ) ")
plt.ylabel("sigma2 = np.logspace(-6, -8, num=20)")
pb.show()



In [241]:

    
vary_baryon = np.linspace(0.0050, 0.050, num=20 ) #set default num = 20
sigma2 = np.logspace(-7, -9, num=20) #set default num = 30

xxx = vary_baryon
yyy = sigma2

zzz = np.array([[LogLikehood_wNoise_1e12_Contour(np.asarray(i), np.asarray(j)) for i in xxx] for j in yyy])

zzzreshaped = zzz.reshape(20,20)

plt.figure()
CS = plt.contour(xxx, yyy, zzzreshaped)
plt.clabel(CS, inline=1, fontsize=10)
plt.xlabel("vary_baryon = np.linspace(0.0050, 0.050, num=20 ) ")
plt.ylabel("sigma2 = np.logspace(-7, -9, num=20)")
pb.show()



In [242]:

    
vary_baryon = np.linspace(0.0050, 0.050, num=20 ) #set default num = 20
sigma2 = np.logspace(-8, -10, num=20) #set default num = 30

xxx = vary_baryon
yyy = sigma2

zzz = np.array([[LogLikehood_wNoise_1e12_Contour(np.asarray(i), np.asarray(j)) for i in xxx] for j in yyy])

zzzreshaped = zzz.reshape(20,20)

plt.figure()
CS = plt.contour(xxx, yyy, zzzreshaped)
plt.clabel(CS, inline=1, fontsize=10)
plt.xlabel("vary_baryon = np.linspace(0.0050, 0.050, num=20 ) ")
plt.ylabel("sigma2 = np.logspace(-8, -10, num=20)")
pb.show()



In [243]:

    
vary_baryon = np.linspace(0.0050, 0.050, num=20 ) #set default num = 20
sigma2 = np.logspace(-9, -11, num=20) #set default num = 30

xxx = vary_baryon
yyy = sigma2

zzz = np.array([[LogLikehood_wNoise_1e12_Contour(np.asarray(i), np.asarray(j)) for i in xxx] for j in yyy])

zzzreshaped = zzz.reshape(20,20)

plt.figure()
CS = plt.contour(xxx, yyy, zzzreshaped)
plt.clabel(CS, inline=1, fontsize=10)
plt.xlabel("vary_baryon = np.linspace(0.0050, 0.050, num=20 ) ")
plt.ylabel("sigma2 = np.logspace(-9, -11, num=20)")
pb.show()



In [244]:

    
vary_baryon = np.linspace(0.0050, 0.050, num=20 ) #set default num = 20
sigma2 = np.logspace(-10, -12, num=20) #set default num = 30

xxx = vary_baryon
yyy = sigma2

zzz = np.array([[LogLikehood_wNoise_1e12_Contour(np.asarray(i), np.asarray(j)) for i in xxx] for j in yyy])

zzzreshaped = zzz.reshape(20,20)

plt.figure()
CS = plt.contour(xxx, yyy, zzzreshaped)
plt.clabel(CS, inline=1, fontsize=10)
plt.xlabel("vary_baryon = np.linspace(0.0050, 0.050, num=20 ) ")
plt.ylabel("sigma2 = np.logspace(-10, -12, num=20)")
pb.show()



In [245]:

    
vary_baryon = np.linspace(0.0050, 0.050, num=20 ) #set default num = 20
sigma2 = np.logspace(-11, -13, num=20) #set default num = 30

xxx = vary_baryon
yyy = sigma2

zzz = np.array([[LogLikehood_wNoise_1e12_Contour(np.asarray(i), np.asarray(j)) for i in xxx] for j in yyy])

zzzreshaped = zzz.reshape(20,20)

plt.figure()
CS = plt.contour(xxx, yyy, zzzreshaped)
plt.clabel(CS, inline=1, fontsize=10)
plt.xlabel("vary_baryon = np.linspace(0.0050, 0.050, num=20 ) ")
plt.ylabel("sigma2 = np.logspace(-11, -13, num=20)")
pb.show()



In [246]:

    
vary_baryon = np.linspace(0.0050, 0.050, num=20 ) #set default num = 20
sigma2 = np.logspace(-12, -14, num=20) #set default num = 30

xxx = vary_baryon
yyy = sigma2

zzz = np.array([[LogLikehood_wNoise_1e12_Contour(np.asarray(i), np.asarray(j)) for i in xxx] for j in yyy])

zzzreshaped = zzz.reshape(20,20)

plt.figure()
CS = plt.contour(xxx, yyy, zzzreshaped)
plt.clabel(CS, inline=1, fontsize=10)
plt.xlabel("vary_baryon = np.linspace(0.0050, 0.050, num=20 ) ")
plt.ylabel("sigma2 = np.logspace(-12, -14, num=20)")
pb.show()



In [247]:

    
vary_baryon = np.linspace(0.0050, 0.050, num=20 ) #set default num = 20
sigma2 = np.logspace(-13, -15, num=20) #set default num = 30

xxx = vary_baryon
yyy = sigma2

zzz = np.array([[LogLikehood_wNoise_1e12_Contour(np.asarray(i), np.asarray(j)) for i in xxx] for j in yyy])

zzzreshaped = zzz.reshape(20,20)

plt.figure()
CS = plt.contour(xxx, yyy, zzzreshaped)
plt.clabel(CS, inline=1, fontsize=10)
plt.xlabel("vary_baryon = np.linspace(0.0050, 0.050, num=20 ) ")
plt.ylabel("sigma2 = np.logspace(-13, -15, num=20)")
pb.show()



In [248]:

    
vary_baryon = np.linspace(0.0050, 0.050, num=20 ) #set default num = 20
sigma2 = np.logspace(-14, -16, num=20) #set default num = 30

xxx = vary_baryon
yyy = sigma2

zzz = np.array([[LogLikehood_wNoise_1e12_Contour(np.asarray(i), np.asarray(j)) for i in xxx] for j in yyy])

zzzreshaped = zzz.reshape(20,20)

plt.figure()
CS = plt.contour(xxx, yyy, zzzreshaped)
plt.clabel(CS, inline=1, fontsize=10)
plt.xlabel("vary_baryon = np.linspace(0.0050, 0.050, num=20 ) ")
plt.ylabel("sigma2 = np.logspace(-14, -16, num=20)")
pb.show()



In [249]:

    
vary_baryon = np.linspace(0.0050, 0.050, num=20) #set default num = 20
sigma2 = np.logspace(-15, -17, num=20) #set default num = 30

xxx = vary_baryon
yyy = sigma2

zzz = np.array([[LogLikehood_wNoise_1e12_Contour(np.asarray(i), np.asarray(j)) for i in xxx] for j in yyy])

zzzreshaped = zzz.reshape(20,20)

plt.figure()
CS = plt.contour(xxx, yyy, zzzreshaped)
plt.clabel(CS, inline=1, fontsize=10)
plt.xlabel("vary_baryon = np.linspace(0.0050, 0.050, num=20 ) ")
plt.ylabel("sigma2 = np.logspace(-15, -17, num=20)")
pb.show()



In [250]:

    
vary_baryon = np.linspace(0.0050, 0.050, num=20 ) #set default num = 20
sigma2 = np.logspace(-16, -18, num=20) #set default num = 30

xxx = vary_baryon
yyy = sigma2

zzz = np.array([[LogLikehood_wNoise_1e12_Contour(np.asarray(i), np.asarray(j)) for i in xxx] for j in yyy])

zzzreshaped = zzz.reshape(20,20)

plt.figure()
CS = plt.contour(xxx, yyy, zzzreshaped)
plt.clabel(CS, inline=1, fontsize=10)
plt.xlabel("vary_baryon = np.linspace(0.0050, 0.050, num=20 ) ")
plt.ylabel("sigma2 = np.logspace(-16, -18, num=20)")
pb.show()



In [251]:

    
vary_baryon = np.linspace(0.010, 0.040, num=30) #set default num = 20
sigma2 = np.logspace(-8, -10, num=30) #set default num = 30

xxx = vary_baryon
yyy = sigma2

zzz = np.array([[LogLikehood_wNoise_1e12_Contour(np.asarray(i), np.asarray(j)) for i in xxx] for j in yyy])

zzzreshaped = zzz.reshape(30,30)

plt.figure()
CS = plt.contour(xxx, yyy, zzzreshaped)
plt.clabel(CS, inline=1, fontsize=10)
plt.xlabel("vary_baryon = np.linspace(0.010, 0.040, num=30)")
plt.ylabel("sigma2 = np.logspace(-8, -10, num=30)")
pb.show()



In [252]:

    
#
# Hold sigma^2 constant, vary C3
#
# Begin at sigma^2 = 1e-5
#
print(vary_baryon)









    



[ 0.01        0.01103448  0.01206897  0.01310345  0.01413793  0.01517241
  0.0162069   0.01724138  0.01827586  0.01931034  0.02034483  0.02137931
  0.02241379  0.02344828  0.02448276  0.02551724  0.02655172  0.02758621
  0.02862069  0.02965517  0.03068966  0.03172414  0.03275862  0.0337931
  0.03482759  0.03586207  0.03689655  0.03793103  0.03896552  0.04      ]



In [253]:

    
tempp = (1e6)*patch # multiply CMB maps by 1e6

def LogLikehood_wNoise_constant(param, sig):
    # param is our parameter, C_3
    Sij = param[:, None, None] * correctmatrix[None, :, :]
    newSij = (1e12)*Sij   # multiply S_ij by 1e12
    Nij = sig * id_mat[None, :, :]
    newNij = (1e12)*Nij
    # Format 7/4pi * param * P_3(M) where param is the parameter we vary, C_l
    # Sij.shape = (20, 3072, 3072)
    Cij = newSij + newNij
    #invCij = np.linalg.inv(Cij)
    logdetC = np.linalg.slogdet(Cij)  # returns sign and determinant; use logdetC[1]
    # model_fit_terms = m^T C^-1 m
    #
    # model_fit_terms = np.array([np.dot(tempval.T , np.dot(invCij[i] , tempval) ) 
    # for i in range(invCij.shape[0])])
    #
    model_fit_terms = np.array([np.dot(tempp.T , np.linalg.solve(Cij[i], tempp) ) for i in range(Cij.shape[0]) ]) 
    return model_fit_terms + logdetC[1] + Npix2pi


vary_baryon = np.linspace(0.0050, 0.040, num=50) #set default num = 20
sigmatwo = 1e-8
print(vary_baryon)

xx = vary_baryon
yy = LogLikehood_wNoise_constant(vary_baryon, sigmatwo)

plt.figure()
CS = plt.plot(xx, yy)
plt.xlabel("vary_baryon = np.linspace(0.010, 0.040, num=30), sigmatwo = 1e-8")
plt.ylabel("LogLF")
pb.show()









    



[ 0.005       0.00571429  0.00642857  0.00714286  0.00785714  0.00857143
  0.00928571  0.01        0.01071429  0.01142857  0.01214286  0.01285714
  0.01357143  0.01428571  0.015       0.01571429  0.01642857  0.01714286
  0.01785714  0.01857143  0.01928571  0.02        0.02071429  0.02142857
  0.02214286  0.02285714  0.02357143  0.02428571  0.025       0.02571429
  0.02642857  0.02714286  0.02785714  0.02857143  0.02928571  0.03
  0.03071429  0.03142857  0.03214286  0.03285714  0.03357143  0.03428571
  0.035       0.03571429  0.03642857  0.03714286  0.03785714  0.03857143
  0.03928571  0.04      ]



In [254]:

    
tempp = (1e6)*patch # multiply CMB maps by 1e6

def LogLikehood_wNoise_constant(param, sig):
    # param is our parameter, C_3
    Sij = param[:, None, None] * correctmatrix[None, :, :]
    newSij = (1e12)*Sij   # multiply S_ij by 1e12
    Nij = sig[:, None, None] * id_mat[None, :, :]
    newNij = (1e12)*Nij
    # Format 7/4pi * param * P_3(M) where param is the parameter we vary, C_l
    # Sij.shape = (20, 3072, 3072)
    Cij = newSij + newNij
    #invCij = np.linalg.inv(Cij)
    logdetC = np.linalg.slogdet(Cij)  # returns sign and determinant; use logdetC[1]
    # model_fit_terms = m^T C^-1 m
    #
    # model_fit_terms = np.array([np.dot(tempval.T , np.dot(invCij[i] , tempval) ) 
    # for i in range(invCij.shape[0])])
    #
    model_fit_terms = np.array([np.dot(tempp.T , np.linalg.solve(Cij[i], tempp) ) for i in range(Cij.shape[0]) ]) 
    return model_fit_terms + logdetC[1] + Npix2pi


vary_baryon = np.linspace(0.0050, 0.040, num=50) #set default num = 20
sigmatwo = np.logspace(-8, -10, num=50)
print(vary_baryon)

xx = vary_baryon
yy = LogLikehood_wNoise_constant(vary_baryon, sigmatwo)

plt.figure()
CS = plt.plot(xx, yy)
plt.xlabel("vary_baryon = np.linspace(0.010, 0.040, num=30), hold sigma2 = np.logspace(-8, -10, num=50)")
plt.ylabel("LogLF")
pb.show()









    



[ 0.005       0.00571429  0.00642857  0.00714286  0.00785714  0.00857143
  0.00928571  0.01        0.01071429  0.01142857  0.01214286  0.01285714
  0.01357143  0.01428571  0.015       0.01571429  0.01642857  0.01714286
  0.01785714  0.01857143  0.01928571  0.02        0.02071429  0.02142857
  0.02214286  0.02285714  0.02357143  0.02428571  0.025       0.02571429
  0.02642857  0.02714286  0.02785714  0.02857143  0.02928571  0.03
  0.03071429  0.03142857  0.03214286  0.03285714  0.03357143  0.03428571
  0.035       0.03571429  0.03642857  0.03714286  0.03785714  0.03857143
  0.03928571  0.04      ]



In [255]:

    
tempp = (1e6)*patch # multiply CMB maps by 1e6

def LogLikehood_wNoise_constant(param, sig):
    # param is our parameter, C_3
    Sij = param[:, None, None] * correctmatrix[None, :, :]
    newSij = (1e12)*Sij   # multiply S_ij by 1e12
    Nij = sig[:, None, None] * id_mat[None, :, :]
    newNij = (1e12)*Nij
    # Format 7/4pi * param * P_3(M) where param is the parameter we vary, C_l
    # Sij.shape = (20, 3072, 3072)
    Cij = newSij + newNij
    #invCij = np.linalg.inv(Cij)
    logdetC = np.linalg.slogdet(Cij)  # returns sign and determinant; use logdetC[1]
    # model_fit_terms = m^T C^-1 m
    #
    # model_fit_terms = np.array([np.dot(tempval.T , np.dot(invCij[i] , tempval) ) 
    # for i in range(invCij.shape[0])])
    #
    model_fit_terms = np.array([np.dot(tempp.T , np.linalg.solve(Cij[i], tempp) ) for i in range(Cij.shape[0]) ]) 
    return model_fit_terms + logdetC[1] + Npix2pi


vary_baryon = np.linspace(0.0050, 0.040, num=50) #set default num = 20
sigmatwo = np.logspace(-7, -10, num=50)
print(vary_baryon)

xx = vary_baryon
yy = LogLikehood_wNoise_constant(vary_baryon, sigmatwo)

plt.figure()
CS = plt.plot(xx, yy)
plt.xlabel("vary_baryon = np.linspace(0.010, 0.040, num=30), sigmatwo = np.logspace(-7, -10, num=50)")
plt.ylabel("LogLF")
pb.show()









    



[ 0.005       0.00571429  0.00642857  0.00714286  0.00785714  0.00857143
  0.00928571  0.01        0.01071429  0.01142857  0.01214286  0.01285714
  0.01357143  0.01428571  0.015       0.01571429  0.01642857  0.01714286
  0.01785714  0.01857143  0.01928571  0.02        0.02071429  0.02142857
  0.02214286  0.02285714  0.02357143  0.02428571  0.025       0.02571429
  0.02642857  0.02714286  0.02785714  0.02857143  0.02928571  0.03
  0.03071429  0.03142857  0.03214286  0.03285714  0.03357143  0.03428571
  0.035       0.03571429  0.03642857  0.03714286  0.03785714  0.03857143
  0.03928571  0.04      ]



In [256]:

    
tempp = (1e6)*patch # multiply CMB maps by 1e6

def LogLikehood_wNoise_constant(param, sig):
    # param is our parameter, C_3
    Sij = param[:, None, None] * correctmatrix[None, :, :]
    newSij = (1e12)*Sij   # multiply S_ij by 1e12
    Nij = sig[:, None, None] * id_mat[None, :, :]
    newNij = (1e12)*Nij
    # Format 7/4pi * param * P_3(M) where param is the parameter we vary, C_l
    # Sij.shape = (20, 3072, 3072)
    Cij = newSij + newNij
    #invCij = np.linalg.inv(Cij)
    logdetC = np.linalg.slogdet(Cij)  # returns sign and determinant; use logdetC[1]
    # model_fit_terms = m^T C^-1 m
    #
    # model_fit_terms = np.array([np.dot(tempval.T , np.dot(invCij[i] , tempval) ) 
    # for i in range(invCij.shape[0])])
    #
    model_fit_terms = np.array([np.dot(tempp.T , np.linalg.solve(Cij[i], tempp) ) for i in range(Cij.shape[0]) ]) 
    return model_fit_terms + logdetC[1] + Npix2pi


vary_baryon = np.linspace(0.0050, 0.040, num=50) #set default num = 20
sigmatwo = np.logspace(-7, -9, num=50)
print(vary_baryon)

xx = vary_baryon
yy = LogLikehood_wNoise_constant(vary_baryon, sigmatwo)

plt.figure()
CS = plt.plot(xx, yy)
plt.xlabel("vary_baryon = np.linspace(0.010, 0.040, num=30), sigmatwo = np.logspace(-7, -9, num=50)")
plt.ylabel("LogLF")
pb.show()









    



[ 0.005       0.00571429  0.00642857  0.00714286  0.00785714  0.00857143
  0.00928571  0.01        0.01071429  0.01142857  0.01214286  0.01285714
  0.01357143  0.01428571  0.015       0.01571429  0.01642857  0.01714286
  0.01785714  0.01857143  0.01928571  0.02        0.02071429  0.02142857
  0.02214286  0.02285714  0.02357143  0.02428571  0.025       0.02571429
  0.02642857  0.02714286  0.02785714  0.02857143  0.02928571  0.03
  0.03071429  0.03142857  0.03214286  0.03285714  0.03357143  0.03428571
  0.035       0.03571429  0.03642857  0.03714286  0.03785714  0.03857143
  0.03928571  0.04      ]



In [257]:

    
vary_baryon = np.linspace(-0.010, 0.040, num=30) #set default num = 20
sigma2 = np.logspace(-8, -10, num=30) #set default num = 30

xxx = vary_baryon
yyy = sigma2

zzz = np.array([[LogLikehood_wNoise_1e12_Contour(np.asarray(i), np.asarray(j)) for i in xxx] for j in yyy])

zzzreshaped = zzz.reshape(30,30)

plt.figure()
CS = plt.contour(xxx, yyy, zzzreshaped)
plt.clabel(CS, inline=1, fontsize=10)
plt.xlabel("vary_baryon = np.linspace(-0.010, 0.040, num=30)")
plt.ylabel("sigma2 = np.logspace(-8, -10, num=30)")
pb.show()



In [ ]:



In [ ]:



In [ ]:



In [ ]:



In [ ]: