

In [1]:

    
%matplotlib inline
import matplotlib.pyplot as plt
import numpy as np
from scipy.io import FortranFile

The following block builds the TT binning matrix (W) following the formula in the Planck 2015 likelihood paper. The Wprime is the binning matrix that binns twice as much (except for the upper end of the scale, there where an odd number of binns.)



In [2]:

    
W=np.zeros([215,2479])
l=30
for b in range(215):
    s=0
    if b < 14:
        for i in range(l,l+5):
            s+=i*(i+1)
        m=l
        for l in range(m,m+5):
            W[b,l-30]=(l*(l+1))/float(s)
        l=l+1
    elif b<170:
        for i in range(l,l+9):
            s+=i*(i+1)
        m=l
        for l in range(m,m+9):
            W[b,l-30]=l*(l+1)/float(s)
        l=l+1
    elif b<200:
        for i in range(l,l+17):
            s+=i*(i+1)
        m=l
        for l in range(m,m+17):
            W[b,l-30]=l*(l+1)/float(s)
        l=l+1
    elif b<215:
        for i in range(l,l+33):
            s+=i*(i+1)
        m=l
        for l in range(m,m+33):
            W[b,l-30]=l*(l+1)/float(s)
        l=l+1
    else:
        print "Fail"
Wprime=np.zeros([115,2479])
l=30
for b in range(115):
    s=0
    if b < 7:
        for i in range(l,l+10):
            s+=i*(i+1)
        m=l
        for l in range(m,m+10):
            Wprime[b,l-30]=l*(l+1)/float(s)
        l=l+1
    elif b<85:
        for i in range(l,l+18):
            s+=i*(i+1)
        m=l
        for l in range(m,m+18):
            Wprime[b,l-30]=l*(l+1)/float(s)
        l=l+1
    elif b<100:
        for i in range(l,l+34):
            s+=i*(i+1)
        m=l
        for l in range(m,m+34):
            Wprime[b,l-30]=l*(l+1)/float(s)
        l=l+1
    elif b<115:
        for i in range(l,l+33):
            s+=i*(i+1)
        m=l
        for l in range(m,m+33):
            Wprime[b,l-30]=l*(l+1)/float(s)
        l=l+1

Load the derivative vectors from file



In [32]:

    
dcls=np.loadtxt("xmat_cosines.txt")

Load, rearrange, and slice the Planck data inverse covariance matrix



In [4]:

    
f=FortranFile('cinv','r')
cinv=f.read_reals(dtype='float').reshape(613,613)



In [5]:

    
for i in range (cinv.shape[0]):
    for j in range(i):
        cinv[j,i]=cinv[i,j]



In [6]:

    
cinv=cinv[0:215,0:215]



In [7]:

    
#need to rebinn the cinv for the prime case
#"unbinning"- use the pseudo inverse of W
rebinn=np.dot(Wprime,np.linalg.pinv(W))
cinvp=np.dot(rebinn,np.dot(cinv,rebinn.T))



In [35]:

    
dclsb=np.dot(W,dcls)
dclsbp=np.dot(Wprime,dcls)



In [36]:

    
evals,evecs=np.linalg.eigh(np.dot(dclsb.T,np.dot(cinv,dclsb)))
evalsp,evecsp=np.linalg.eigh(np.dot(dclsbp.T,np.dot(cinvp,dclsbp)))



In [40]:

    
plt.gca().set_yscale('log')
plt.plot(evals,".")
plt.show()



In [41]:

    
plt.gca().set_yscale('log')
plt.plot(evalsp,".")
plt.show()



In [52]:

    
plt.plot(evecs[:,4])
plt.show()



In [53]:

    
plt.plot(evecsp[:,4])
plt.show()



In [67]:

    
plt.plot(evecs[:,4]-evecsp[:,4])
plt.show()



In [60]:

    
plt.plot(np.dot(dcls,evecs[:,4].T))
plt.show()



In [59]:

    
plt.plot(np.dot(dcls,evecsp[:,4].T))
plt.show()



In [66]:

    
plt.plot(np.dot(dcls,(evecs[:,4].T-evecsp[:,4].T)))
plt.show()



In [ ]: