In [1]:

    

#
# LF test with noise
# C = S + N, where noise N = sigma^2 I, where I is the identity matrix
#
# We generate 20 values of parameters x for function logLF(x)
# We chose our parameters around the value of Boltzmann code generated C_3
# CAMB generated C_3 equals 5.88275e-10
#
# 1. Set all C_l = 0 except C_3  DONE
#
# 2. LF is based on a matrix C that only has a P_3 term. DONE
#
# 3. Add to that matrix a white noise term (sigma^2 on the diagonal). DONE
#
# 4. LF now has *exactly two* free parameters, C_3 and sigma^2. DONE
#
# 5. What is LF vs C_3 (at, perhaps, a couple values of sigma^2)? 
#
# 6. over-plot |a_3m|^2 values
#
#


    

In [2]:

    

%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 [3]:

    

#
# Review of Likelihood Formalism:
#
# -2 loglikelihood is
# -2 ln L \propto m^T C^-1 m + lnDet C
# where C = S + N
# We are working with noiseless maps, N = 0, so C = S
#
# In real space, 
# data: the temperature map, 
# parameters: theoretical CAMB generated C_l, C^theory_l 
#
# m = array of temperature pixels
# S = S_ij
# N = diagonal noise (but here it is noiseless! N=0)
#
# In spherical harmonic space
# data: the Healpix anafast outputarray of a_lm (or similarly, \hat{C}_l )
# parameters: the theoretical C_l
#
# m = a_lm coefficients
# S is diagonal (C_2, C_3, etc.)
# N is non-sparse matrix (but here it is noiseless! N=0)
#
#
# NOTE: The correct standard is to systematically remove the monopole, dipole terms l=0,l=1
#       Also, we use in the following lmax = 2*nside


    

In [4]:

    

cd ~/Desktop/CMBintheLikeHoodz/Likelihood_Comparison


    

    




/Users/evanbiederstedt/Desktop/CMBintheLikeHoodz/Likelihood_Comparison

In [5]:

    

camb1 = "camb_nside16_lmax32_alms.fits"
camb2 = "camb_nside16_lmax32_map.fits"
camb3 = "camb_nside16_lmax32_scalcls.fits"
planck1 = "100GHz_nside16_lmax32_cls.fits"
planck2 = "100GHz_nside16_lmax32_cmb_alm.fits"
planck3 = "100GHz_nside16_lmax32_sky_alm.fits"
planck4 = "100GHz_nside16_lmax32_skymap.fits"

nside = 16


    

In [6]:

    

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 3072
The maximum ell of the power spectrum C_l set to lmax = 2*nside 32
Healpix tells me total number of pixels npix is equal to 3072

In [7]:

    

#
# Begin with a Munich Planck-simulated map, and CAMB Boltzmann-code generated C_l values
#


    

In [8]:

    

# Theoretical scalar C_l array, CAMB
#
# open a FITS file, theoretical C_l values generated by CAMB
# type()=pyfits.hdu.hdulist.HDUList
cl_open = pf.open(camb3)
# recall camb3 = "camb_nside16_lmax32_scalcls.fits"


    

In [9]:

    

theoryCls_arr1 = cl_open[1].data
print theoryCls_arr1[:10]
# Recall there are four columns: temp, E pol, B pol, grad-temp cross terms
# first two values are zero, i.e. monopole, dipole
# XXX.field() references columns by 0-index
# field(0) is temperature values
# all Cl scalar temp values put into ndarray
# type()=numpy.ndarray


    

    




[(0.0, 0.0, 0.0, 0.0) (0.0, 0.0, 0.0, 0.0)
 (1.2639207e-09, 5.5308927e-14, 0.0, 3.6016001e-12)
 (5.8827504e-10, 4.5345086e-14, 0.0, 2.2704995e-12)
 (3.2867314e-10, 3.0680604e-14, 0.0, 1.4030458e-12)
 (2.065753e-10, 1.8198442e-14, 0.0, 8.6680196e-13)
 (1.4100289e-10, 9.6222612e-15, 0.0, 5.4035639e-13)
 (1.0233801e-10, 4.5884262e-15, 0.0, 3.3960166e-13)
 (7.7812937e-11, 2.0448748e-15, 0.0, 2.1437516e-13)
 (6.1362422e-11, 9.4503657e-16, 0.0, 1.3715881e-13)]

In [10]:

    

cltemps = theoryCls_arr1.field(0) 
print cltemps
print "The length of the array of theoretical Cl's is " +str(len(cltemps))
print "The array contains [C_0, C_1, C_2,..., C_" +str(len(cltemps)-1) + "]"
#print type(cltemps)=numpy.ndarray


    

    




[  0.00000000e+00   0.00000000e+00   1.26392075e-09   5.88275040e-10
   3.28673144e-10   2.06575299e-10   1.41002890e-10   1.02338013e-10
   7.78129366e-11   6.13624221e-11   4.98225766e-11   4.14123076e-11
   3.50811741e-11   3.01814025e-11   2.63133803e-11   2.32104266e-11
   2.06850804e-11   1.86004823e-11   1.68565179e-11   1.53792031e-11
   1.41130128e-11   1.30162998e-11   1.20591913e-11   1.12186076e-11
   1.04760627e-11   9.81660760e-12   9.22804801e-12   8.70035207e-12
   8.22520958e-12   7.79568684e-12   7.40596039e-12   7.05096528e-12
   6.72585642e-12]
The length of the array of theoretical Cl's is 33
The array contains [C_0, C_1, C_2,..., C_32]

In [11]:

    

# remove monopole l=0 and dipole l=1
theoryCl = cltemps[2:]
# len(theoryCl) = 31
print theoryCl
# theoryCl is np.ndarray of theoretical [C_2, C_3, C_4, ..., C_32]


    

    




[  1.26392075e-09   5.88275040e-10   3.28673144e-10   2.06575299e-10
   1.41002890e-10   1.02338013e-10   7.78129366e-11   6.13624221e-11
   4.98225766e-11   4.14123076e-11   3.50811741e-11   3.01814025e-11
   2.63133803e-11   2.32104266e-11   2.06850804e-11   1.86004823e-11
   1.68565179e-11   1.53792031e-11   1.41130128e-11   1.30162998e-11
   1.20591913e-11   1.12186076e-11   1.04760627e-11   9.81660760e-12
   9.22804801e-12   8.70035207e-12   8.22520958e-12   7.79568684e-12
   7.40596039e-12   7.05096528e-12   6.72585642e-12]

In [12]:

    

# Our input data is Gaerching generated, noiseless full-sky map
# Temperature map: here we use Planck simulated map from Munich, not CAMB map
# http://gavo.mpa-garching.mpg.de/planck/
#
# Read in with Healpy routine/function
#
# Use planck4 = "100GHz_nside16_lmax32_skymap.fits"
# This is a simulated data, 100GHz (where CMB dominates), no foregrounds
# 
mapread_camb2 = hp.read_map(camb2) # Healpix routine, input the sky map


    

    




NSIDE = 16
ORDERING = RING in fits file

In [13]:

    

hp.mollview(mapread_camb2) # visualization of full-sky CMB map, nside=16, lmax=32


    

In [14]:

    

# The uploaded temperature map is mapread_planck4 = hp.read_map(planck4)
print type(mapread_camb2) # type(mapread_planck4) = np.ndarray
print mapread_camb2.shape # mapread_planck4.shape = (3072, ) = (N_pix, )
#
# rename array for convenience
tempval = mapread_camb2
#print tempval


    

    




<type 'numpy.ndarray'>
(3072,)

In [15]:

    

# Next, we use healpy map2alm to tranform to alm values
# Our input data is CAMB generated, noiseless full-sky map
# We calculate an array of a_lm from this by using Healpix map2alm, a subroutine of anafast
#
# map2alm only outputs m >=0 values, because m = -l values are equivalent to m = +l values
#
# Using map2alm, the length of the alm array is expected to be:
# (mmax * (2 * lmax + 1 - mmax)) / 2 + lmax + 1)"
#
# For mmax = lmax, this is l(l+1)/2 + l + 1
# i.e. 
# l = 0, there is 1
# l = 1, there is 3
# l = 2, there is 6
# l = 3, there is 10
# l = 4, there is 15
# etc.
almarr = hp.map2alm(mapread_camb2)                # This is an array of a_lm values
print "The array of spherical harmonic coefficients a_lm is"
#print almarr
print "The arr.shape is " + str(almarr.shape)
print "The length of a_lm array is " + str(len(almarr))
#
print "For l=3, map2alm gives (a_00, a_10, a_11, a_20, a_21, a_22, a_30, a_31, a_32, a_33)"
print "However, this is NOT the order of the output! See below"
# In the Fortran F90 subroutines, complex alm are stored in an array that has 
# two dimensions to contain coefficients for positive and negative m values.
# Healpy doesn't do this....I think

print "============================="
print "============================="
print "Check indices with healpy.sphtfunc.Alm.getidx(lmax, l, m)"
print "Default ordering of healpy.map2alm() output is "
print "(0,0), (1,0), ..., (lmax, 0),"
print "(1,1), (2,1), ...., (lmax, 1),"
print "(2,2), .... (lmax, 2),(3,3), ...., (lmax, 3), etc. , .... (lmax, lmax)."


    

    




The array of spherical harmonic coefficients a_lm is
The arr.shape is (1176,)
The length of a_lm array is 1176
For l=3, map2alm gives (a_00, a_10, a_11, a_20, a_21, a_22, a_30, a_31, a_32, a_33)
However, this is NOT the order of the output! See below
=============================
=============================
Check indices with healpy.sphtfunc.Alm.getidx(lmax, l, m)
Default ordering of healpy.map2alm() output is 
(0,0), (1,0), ..., (lmax, 0),
(1,1), (2,1), ...., (lmax, 1),
(2,2), .... (lmax, 2),(3,3), ...., (lmax, 3), etc. , .... (lmax, lmax).

In [16]:

    

# ==========================
# DEMONSTRATION
# Notice how a_lm is indexed
# ==========================

mmm = np.arange(12) # define a map, i.e. an array of 12 "pixels"
lmaxxx = 4
alm = hp.map2alm(mmm, lmax=lmaxxx) # spherical harmonic transform
lm = hp.map2alm(mmm, lmax=lmaxxx) # spherical harmonic transform
print(alm)
print(alm.shape)

# So alm is actually a 1D vector.
# How is alm indexed?


l, m = hp.Alm.getlm(lmax=lmaxxx)
print(l)

print(m)




print "The l values are "+str(l)
print "The m values are "+str(m)
print " (l,m) is in order " +str(list(zip(l,m)))
#
# l, m = hp.Alm.getlm(lmax=lmax)
# print(l)
# [0 1 2 1 2 2]
# print(m)
# [0 0 0 1 1 2]
#
#
# So, for l = 2, m is [0, 1, 2].
#
# ==========================
# Notice how a_lm is indexed
# ==========================
#
#
#


    

    




[  1.43186913e+01 +0.00000000e+00j  -9.06158785e+00 +0.00000000e+00j
  -7.57672587e+00 +0.00000000e+00j   5.38292381e+00 +0.00000000e+00j
   2.49127749e+00 +0.00000000e+00j   6.02204636e-01 -2.01587603e+00j
   4.55476051e-16 +1.48937220e-15j   2.43950547e-01 -3.28507898e-01j
  -2.47937026e-17 +3.49331635e-16j   3.65595455e-01 -7.92402132e+00j
   5.92284905e-17 -4.53852723e-15j  -3.16614952e-01 -1.44872968e+01j
  -2.28780515e-03 +4.33449165e-01j  -5.66350430e-16 +1.59192782e-15j
   1.08501635e+01 -3.35394518e-16j]
(15,)
[0 1 2 3 4 1 2 3 4 2 3 4 3 4 4]
[0 0 0 0 0 1 1 1 1 2 2 2 3 3 4]
The l values are [0 1 2 3 4 1 2 3 4 2 3 4 3 4 4]
The m values are [0 0 0 0 0 1 1 1 1 2 2 2 3 3 4]
 (l,m) is in order [(0, 0), (1, 0), (2, 0), (3, 0), (4, 0), (1, 1), (2, 1), (3, 1), (4, 1), (2, 2), (3, 2), (4, 2), (3, 3), (4, 3), (4, 4)]

In [17]:

    

# Check with healpy.sphtfunc.Alm.getidx(lmax, l, m)
# Returns index corresponding to (l,m) in an array describing alm up to lmax.
#
ell, emm = hp.Alm.getlm(lmax=32)

print "len(ell) is " +str(len(ell))
print "len(emm) is "+str(len(emm))
print "l values are "+str(ell[:10])
print "m values are "+str(emm[:10])
pairs = list(zip(ell, emm)) # put values together in pairs, zip()
ellemm = np.vstack((ell,emm)).T  # equivalent to list(zip(ell,emm)), but uses numpy throughout
print "Indices for a_lm for lmax (l, m) are:" 
print str(pairs[:50]) # The expected output


    

    




len(ell) is 561
len(emm) is 561
l values are [0 1 2 3 4 5 6 7 8 9]
m values are [0 0 0 0 0 0 0 0 0 0]
Indices for a_lm for lmax (l, m) are:
[(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), (13, 0), (14, 0), (15, 0), (16, 0), (17, 0), (18, 0), (19, 0), (20, 0), (21, 0), (22, 0), (23, 0), (24, 0), (25, 0), (26, 0), (27, 0), (28, 0), (29, 0), (30, 0), (31, 0), (32, 0), (1, 1), (2, 1), (3, 1), (4, 1), (5, 1), (6, 1), (7, 1), (8, 1), (9, 1), (10, 1), (11, 1), (12, 1), (13, 1), (14, 1), (15, 1), (16, 1), (17, 1)]

In [18]:

    

print ellemm[:10]


    

    




[[0 0]
 [1 0]
 [2 0]
 [3 0]
 [4 0]
 [5 0]
 [6 0]
 [7 0]
 [8 0]
 [9 0]]

In [19]:

    

#
# For our first test, mode l = 3, we need to access a_lm coefficients a_30, a_31, a_32, a_33
# To find this for lmax = 32, we use
# healpy.sphtfunc.Alm.getidx(lmax, l, m)
# Returns index corresponding to (l,m) in an array describing alm up to lmax.
#
# Find the indices
index_a30 = hp.Alm.getidx(lmax=32, l=3, m=0)
index_a31 = hp.Alm.getidx(lmax=32, l=3, m=1)
index_a32 = hp.Alm.getidx(lmax=32, l=3, m=2)
index_a33 = hp.Alm.getidx(lmax=32, l=3, m=3)


    

In [20]:

    

print "Index a_30 is " +str(index_a30)
print "Index a_31 is "+str(index_a31)
print "Index a_32 is "+str(index_a32)
print "Index a_33 is "+str(index_a33)


    

    




Index a_30 is 3
Index a_31 is 35
Index a_32 is 66
Index a_33 is 96

In [21]:

    

#
# Create an array with only the values a_3m, i.e. a_30, a_31, a_32, a_33
#
# First convert the array of alm coefficients into a real
#
realalm = almarr.real
#
print realalm[:36]


    

    




[  4.27365741e-10   2.29584030e-09  -3.99058124e-05  -3.54786292e-05
   2.10319852e-05   9.87750024e-06   7.97092193e-06  -1.18858780e-05
  -7.27840112e-06   3.61680988e-07  -4.90787138e-06   3.47579016e-06
  -4.30349062e-06  -2.65684938e-07  -1.00274382e-06  -1.12408678e-06
   4.89651663e-06   7.38091856e-07   2.10038079e-07   3.19967529e-06
  -2.82808712e-06   3.98389365e-07   9.92859235e-07   1.12017252e-07
   1.85264741e-06  -1.88066911e-06  -5.47687726e-07  -4.03524074e-08
   6.09934135e-07   6.46948184e-07  -2.27504639e-07   6.00326431e-07
  -1.16138884e-07   6.93299985e-09   3.15959247e-09   6.92273981e-09]

In [22]:

    

empty_almlist = []
#
a30 = realalm[3]
a31 = realalm[35]
a32 = realalm[66]
a33 = realalm[96]
#
print "a30 is " + str(a30)
print "a31 is " + str(a31)
print "a32 is " + str(a32)
print "a33 is " + str(a33)
#
print str(pairs[3]) # Check with our output above
print str(pairs[35])
print str(pairs[66])
print str(pairs[96])
#
empty_almlist.append(a30)
empty_almlist.append(a31)
empty_almlist.append(a32)
empty_almlist.append(a33)
#
print empty_almlist


    

    




a30 is -3.54786291692e-05
a31 is 6.92273981371e-09
a32 is -7.65807282676e-07
a33 is 2.69713402623e-05
(3, 0)
(3, 1)
(3, 2)
(3, 3)
[-3.5478629169177805e-05, 6.9227398137090744e-09, -7.6580728267559654e-07, 2.6971340262298287e-05]

In [23]:

    

# create array of real-valued alm coefficients, a30 a31 a32 a33
realalm3 = np.asarray(empty_almlist) # np.asarray() converts input into an array
print realalm3


    

    




[ -3.54786292e-05   6.92273981e-09  -7.65807283e-07   2.69713403e-05]

In [24]:

    

# Repeat the above procedure for mode l = 4, i.e. a40 a41 a42 a43 a44
# Find the indices
index_a40 = hp.Alm.getidx(lmax=32, l=4, m=0)
index_a41 = hp.Alm.getidx(lmax=32, l=4, m=1)
index_a42 = hp.Alm.getidx(lmax=32, l=4, m=2)
index_a43 = hp.Alm.getidx(lmax=32, l=4, m=3)
index_a44 = hp.Alm.getidx(lmax=32, l=4, m=4)
#
print "Index a_40 is " +str(index_a40)
print "Index a_41 is "+str(index_a41)
print "Index a_42 is "+str(index_a42)
print "Index a_43 is "+str(index_a43)
print "Index a_44 is "+str(index_a44)
#
# Check with the above ouput
print str(pairs[4]) 
print str(pairs[36])
print str(pairs[67])
print str(pairs[97])
print str(pairs[126])
#
emptylistalm2 = []
#
#print realalm
#
a40 = realalm[4]
a41 = realalm[36]
a42 = realalm[67]
a43 = realalm[97]
a44 = realalm[127]
#
print "a40 is " + str(a40)
print "a41 is " + str(a41)
print "a42 is " + str(a42)
print "a43 is " + str(a43)
print "a44 is " + str(a44)
#
emptylistalm2.append(a40)
emptylistalm2.append(a41)
emptylistalm2.append(a42)
emptylistalm2.append(a43)
emptylistalm2.append(a44)
#
#print emptylistalm2


    

    




Index a_40 is 4
Index a_41 is 36
Index a_42 is 67
Index a_43 is 97
Index a_44 is 126
(4, 0)
(4, 1)
(4, 2)
(4, 3)
(4, 4)
a40 is 2.10319852146e-05
a41 is 3.20500644244e-09
a42 is 3.94936640444e-07
a43 is 3.29951851271e-06
a44 is -1.19861955538e-09

In [25]:

    

# create array of real-valued alm coefficients, a40 a41 a42 a43 a44
realalm4 = np.asarray(emptylistalm2) # np.asarray() converts input into an array
print realalm4


    

    




[  2.10319852e-05   3.20500644e-09   3.94936640e-07   3.29951851e-06
  -1.19861956e-09]

In [26]:

    

# Calculate (abs(alm))**2 i.e. |alm|^2
abs_alm3 = np.absolute(realalm3)
abs_alm4 = np.absolute(realalm4)
print abs_alm3
print abs_alm4
# Now calculate the squares element-wise, x**2
alm3_squared = abs_alm3**2
alm4_squared = abs_alm4**2
print alm3_squared
print alm4_squared


    

    




[  3.54786292e-05   6.92273981e-09   7.65807283e-07   2.69713403e-05]
[  2.10319852e-05   3.20500644e-09   3.94936640e-07   3.29951851e-06
   1.19861956e-09]
[  1.25873313e-09   4.79243265e-17   5.86460794e-13   7.27453196e-10]
[  4.42344402e-10   1.02720663e-17   1.55974950e-13   1.08868224e-11
   1.43668884e-18]

In [27]:

    

# For l = 3 test, we need theoretical value of C_3; ditto for l = 4
print theoryCl
C3 = theoryCl[1]
print "theory C_3 is " +str(C3)
C4 = theoryCl[2]
print "theory C_4 is "+str(C4)


    

    




[  1.26392075e-09   5.88275040e-10   3.28673144e-10   2.06575299e-10
   1.41002890e-10   1.02338013e-10   7.78129366e-11   6.13624221e-11
   4.98225766e-11   4.14123076e-11   3.50811741e-11   3.01814025e-11
   2.63133803e-11   2.32104266e-11   2.06850804e-11   1.86004823e-11
   1.68565179e-11   1.53792031e-11   1.41130128e-11   1.30162998e-11
   1.20591913e-11   1.12186076e-11   1.04760627e-11   9.81660760e-12
   9.22804801e-12   8.70035207e-12   8.22520958e-12   7.79568684e-12
   7.40596039e-12   7.05096528e-12   6.72585642e-12]
theory C_3 is 5.88275e-10
theory C_4 is 3.28673e-10

In [28]:

    

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


    

    




[ 0  1  2  3  4  5  6  7  8  9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24
 25 26 27 28 29 30 31 32]
[ 2  3  4  5  6  7  8  9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26
 27 28 29 30 31 32]

In [29]:

    

# Calculate an array of (2*l + 1)C_l 
# i.e. 5*C_2, 7*C_3, 9*C_4, 11*C_5, 13*C_6, ...
print theoryCl

for i in ellval:
    paramsCl = (2*ellval + 1)*theoryCl # define array (2*l + 1)C_l 
    
print paramsCl


    

    




[  1.26392075e-09   5.88275040e-10   3.28673144e-10   2.06575299e-10
   1.41002890e-10   1.02338013e-10   7.78129366e-11   6.13624221e-11
   4.98225766e-11   4.14123076e-11   3.50811741e-11   3.01814025e-11
   2.63133803e-11   2.32104266e-11   2.06850804e-11   1.86004823e-11
   1.68565179e-11   1.53792031e-11   1.41130128e-11   1.30162998e-11
   1.20591913e-11   1.12186076e-11   1.04760627e-11   9.81660760e-12
   9.22804801e-12   8.70035207e-12   8.22520958e-12   7.79568684e-12
   7.40596039e-12   7.05096528e-12   6.72585642e-12]
[  6.31960373e-09   4.11792528e-09   2.95805830e-09   2.27232828e-09
   1.83303757e-09   1.53507020e-09   1.32281992e-09   1.16588602e-09
   1.04627411e-09   9.52483075e-10   8.77029352e-10   8.14897868e-10
   7.63088029e-10   7.19523226e-10   6.82607653e-10   6.51016879e-10
   6.23691163e-10   5.99788922e-10   5.78633527e-10   5.59700893e-10
   5.42663609e-10   5.27274559e-10   5.13327074e-10   5.00646987e-10
   4.89086545e-10   4.78519364e-10   4.68836946e-10   4.59945524e-10
   4.51763584e-10   4.44210813e-10   4.37180667e-10]

In [30]:

    

norm = ((2*ellval + 1))/(4*math.pi)
print norm


    

    




[ 0.39788736  0.5570423   0.71619724  0.87535219  1.03450713  1.19366207
  1.35281702  1.51197196  1.6711269   1.83028185  1.98943679  2.14859173
  2.30774667  2.46690162  2.62605656  2.7852115   2.94436645  3.10352139
  3.26267633  3.42183128  3.58098622  3.74014116  3.89929611  4.05845105
  4.21760599  4.37676094  4.53591588  4.69507082  4.85422576  5.01338071
  5.17253565]

In [31]:

    

anafastCl = hp.anafast(mapread_camb2, lmax=32)
#len(anafastCl) = 33
# remove monopole and dipole values, l=0, l=1
hatCl = anafastCl[2:] #len() = 31, type() = np.ndarray
hatC3 = hatCl[1]   # index 0 = C2, 1 = C3, etc.
hatC4 = hatCl[2]
print hatC3
print hatC4


    

    




6.49352544114e-10
3.16267210302e-10

In [32]:

    

# 
# Add a_lm squared, |a_lm|^2
#
print "The values for |a_lm|^2 are : "
print "For |a_3m|**2 such that a_30, a_31, a_32, a_33: "
print str(alm3_squared)
print "And for |a_4m|**2 such that a_40, a_41, a_42, a_43, a_44: "
print str(alm4_squared)


    

    




The values for |a_lm|^2 are : 
For |a_3m|**2 such that a_30, a_31, a_32, a_33: 
[  1.25873313e-09   4.79243265e-17   5.86460794e-13   7.27453196e-10]
And for |a_4m|**2 such that a_40, a_41, a_42, a_43, a_44: 
[  4.42344402e-10   1.02720663e-17   1.55974950e-13   1.08868224e-11
   1.43668884e-18]

In [33]:

    

# =========================================================================
#
# =========================================================================
#
# Data:
#     tempval      # the array of pixel values, (3072,)
#     realalm3     # array of alm values, a30, a31, a32, a33 
#     realalm4     # array of alm values, a40, a41, a42, a43, a44
#     alm3_squared # array of |alm|^2, (abs(a3m))**2
#     alm4_squared # array of |alm|^2, (abs(a4m))**2
#     hatCl        # array of anafast-calculated \hat{C}_l values, l=2 to l=32
#     hatC3        # \hat{C}_3 value
#     hatC4        # \hat{C}_4 value
#
# Parameters:
#     theoryCl     # array of Boltzmann code generated C_l, i.e. C^{theory}_l
#     paramsCl     # array of (2*l + 1)C_l from l=2 to l=lmax
#     C3           # array of C_3 value
#     C4           # array of C_4 value
#    
# Array of ell's: 
#     ellval       # array of l = 2 to l=lmax
#                  # [2 3 4 ... 31 32]
#     norm         # array of (2*l+1)/4pi
#                  # [5/4pi 7/4pi 9/4pi 11/4pi ... 63/4pi 65/4pi]
# =========================================================================
#
# =========================================================================


    

In [34]:

    

#
# 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"
#
# npix = 3072


    

In [35]:

    

totalpix = np.arange(npix) # An array indexing the total number of pixels


    

In [36]:

    

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

    

vecvalx = vecval[0] #shape (3072,)
vecvaly = vecval[1]
vecvalz = vecval[2]


    

In [38]:

    

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

    

trans = totalvecval.T #transpose


    

In [40]:

    

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


    

In [41]:

    

# =========================================================
# =========================================================
#
# \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 [42]:

    

# =========================================================
# =========================================================
#
# Now calculate the 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
#
# STEP 1: do inverse of Sij
#        invSij = np.linalg.inv(Sij)
#
# STEP 2: do matrix mulplication, m.T*inv(C)*m
#       
# lnL_modelfit_terms = np.array([np.dot(tempval.T , np.dot(invSij[i] , tempval) ) for i in range(invSij.shape[0])])
#
# STEP 3: do logdet Sij
#
# logdetC = np.linalg.slogdet(Sij) #computes sign and log det C
# logdetC[1]
#
# STEP 4: compute N_pix * 2pi
#
# Npix2pi = (npix)*2*math.pi
#
# Step 5: -2loglikelihood = m.T*inv(C)*m + logdetC[1] + Npix2pi
# 
# =========================================================
# =========================================================


    

In [43]:

    

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

    

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

    

# JULY 1, 2015
#
#
# Here we define the LF 
#
# Next we should writen the REAL-space log-likelihood, -2 lnLikefunction
# this is a function of parameters
# plot against parameters, theory C_l
#
# Likelihood function should use -2lnL /propto m^T C^-1 M + log det M
#


    

In [46]:

    

# TEST C_3 Likelihood
#
# The covariance matrix is a function of variable "x" where "x" is "C_3", an unknown parameter. 
# 
# Our covariance matrix is therefore S_ij  = 7/4pi * x * P_3(matrix)   
# (We set l=3, i.e. l=3, and P_3)
# 
# The LF is then a function of x, LF(x). 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 [47]:

    

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


    

In [48]:

    

print C3 # Boltzmann code CAMB output


    

    




5.88275e-10

In [49]:

    

# generate a number of samples to plot for x
# Boltzmann code value for C3 is 5.88275e-10
# start at 1e-10, end at 9e-10
# default samples generated is 50


    

In [50]:

    

vary_x_samples1 = np.linspace(1e-10, 9e-10, num=20 ) #set default num = 20


    

In [51]:

    

print vary_x_samples1


    

    




[  1.00000000e-10   1.42105263e-10   1.84210526e-10   2.26315789e-10
   2.68421053e-10   3.10526316e-10   3.52631579e-10   3.94736842e-10
   4.36842105e-10   4.78947368e-10   5.21052632e-10   5.63157895e-10
   6.05263158e-10   6.47368421e-10   6.89473684e-10   7.31578947e-10
   7.73684211e-10   8.15789474e-10   8.57894737e-10   9.00000000e-10]

In [52]:

    

# create Sij array, len()=31, l=2 to l=32
# Sij = norm_matrix * theoryCl[:, None, None]
# [5/4pi*C_2*P_2(M)  7/4pi*C_3*P_3(M) ....   65/4pi*C_32*P_32(M)]


    

In [53]:

    

print hatC3  # anafast generated C_l for l = 3 from sky map
print C3     # CAMB generated C_3, Boltzmann code output


    

    




6.49352544114e-10
5.88275e-10

In [54]:

    

print alm3_squared  # a_3m extracted from sky map, absolute value, squared


    

    




[  1.25873313e-09   4.79243265e-17   5.86460794e-13   7.27453196e-10]

In [55]:

    

np.set_printoptions(threshold=100000)  # Default is threshold=1000
## Use this to print all values, disables corner printing


    

In [56]:

    

# print out all matrices and parts of the logLF
print vary_x_samples1  # generated from 1e-10 to 9e-10


    

    




[  1.00000000e-10   1.42105263e-10   1.84210526e-10   2.26315789e-10
   2.68421053e-10   3.10526316e-10   3.52631579e-10   3.94736842e-10
   4.36842105e-10   4.78947368e-10   5.21052632e-10   5.63157895e-10
   6.05263158e-10   6.47368421e-10   6.89473684e-10   7.31578947e-10
   7.73684211e-10   8.15789474e-10   8.57894737e-10   9.00000000e-10]

In [57]:

    

print Npix2pi # LogLF constant


    

    




19301.9452637

In [58]:

    

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


    

    




[  1.00000000e-12   8.53167852e-13   7.27895384e-13   6.21016942e-13
   5.29831691e-13   4.52035366e-13   3.85662042e-13   3.29034456e-13
   2.80721620e-13   2.39502662e-13   2.04335972e-13   1.74332882e-13
   1.48735211e-13   1.26896100e-13   1.08263673e-13   9.23670857e-14
   7.88046282e-14   6.72335754e-14   5.73615251e-14   4.89390092e-14
   4.17531894e-14   3.56224789e-14   3.03919538e-14   2.59294380e-14
   2.21221629e-14   1.88739182e-14   1.61026203e-14   1.37382380e-14
   1.17210230e-14   1.00000000e-14]

In [59]:

    

print 7*hatC3
print hatC3


    

    




4.5454678088e-09
6.49352544114e-10

In [60]:

    

# Second test with rescaled values
vary_x_samples3 = np.logspace(-8, -12, num=30) #num = 30 
sigma3 = np.logspace(-12, -14, num=30)


    

In [61]:

    

print vary_x_samples3
print sigma3


    

    




[  1.00000000e-08   7.27895384e-09   5.29831691e-09   3.85662042e-09
   2.80721620e-09   2.04335972e-09   1.48735211e-09   1.08263673e-09
   7.88046282e-10   5.73615251e-10   4.17531894e-10   3.03919538e-10
   2.21221629e-10   1.61026203e-10   1.17210230e-10   8.53167852e-11
   6.21016942e-11   4.52035366e-11   3.29034456e-11   2.39502662e-11
   1.74332882e-11   1.26896100e-11   9.23670857e-12   6.72335754e-12
   4.89390092e-12   3.56224789e-12   2.59294380e-12   1.88739182e-12
   1.37382380e-12   1.00000000e-12]
[  1.00000000e-12   8.53167852e-13   7.27895384e-13   6.21016942e-13
   5.29831691e-13   4.52035366e-13   3.85662042e-13   3.29034456e-13
   2.80721620e-13   2.39502662e-13   2.04335972e-13   1.74332882e-13
   1.48735211e-13   1.26896100e-13   1.08263673e-13   9.23670857e-14
   7.88046282e-14   6.72335754e-14   5.73615251e-14   4.89390092e-14
   4.17531894e-14   3.56224789e-14   3.03919538e-14   2.59294380e-14
   2.21221629e-14   1.88739182e-14   1.61026203e-14   1.37382380e-14
   1.17210230e-14   1.00000000e-14]

In [62]:

    

# For N matrix, set the identity
id_mat = np.identity(3072)
# print id_mat    # This is a (3072, 3072) matrix


    

In [63]:

    

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


    

In [64]:

    

# Now, rescale
#
# Multiply the CMB maps by 1e6 and the C matrix entries by 1e12 and try again.
#
tempp = (1e6)*tempval # multiply CMB maps by 1e6

def LogLikehood_wNoise_1e12(param, sig):
    # param is our parameter, C_3
    Sij = param[:, None, None] * norm_matrix[1][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 [ ]:

    




    

In [65]:

    

print sigma3
print "*************"
print vary_x_samples3


    

    




[  1.00000000e-12   8.53167852e-13   7.27895384e-13   6.21016942e-13
   5.29831691e-13   4.52035366e-13   3.85662042e-13   3.29034456e-13
   2.80721620e-13   2.39502662e-13   2.04335972e-13   1.74332882e-13
   1.48735211e-13   1.26896100e-13   1.08263673e-13   9.23670857e-14
   7.88046282e-14   6.72335754e-14   5.73615251e-14   4.89390092e-14
   4.17531894e-14   3.56224789e-14   3.03919538e-14   2.59294380e-14
   2.21221629e-14   1.88739182e-14   1.61026203e-14   1.37382380e-14
   1.17210230e-14   1.00000000e-14]
*************
[  1.00000000e-08   7.27895384e-09   5.29831691e-09   3.85662042e-09
   2.80721620e-09   2.04335972e-09   1.48735211e-09   1.08263673e-09
   7.88046282e-10   5.73615251e-10   4.17531894e-10   3.03919538e-10
   2.21221629e-10   1.61026203e-10   1.17210230e-10   8.53167852e-11
   6.21016942e-11   4.52035366e-11   3.29034456e-11   2.39502662e-11
   1.74332882e-11   1.26896100e-11   9.23670857e-12   6.72335754e-12
   4.89390092e-12   3.56224789e-12   2.59294380e-12   1.88739182e-12
   1.37382380e-12   1.00000000e-12]

In [66]:

    

print vary_x_samples3
print sigma3


    

    




[  1.00000000e-08   7.27895384e-09   5.29831691e-09   3.85662042e-09
   2.80721620e-09   2.04335972e-09   1.48735211e-09   1.08263673e-09
   7.88046282e-10   5.73615251e-10   4.17531894e-10   3.03919538e-10
   2.21221629e-10   1.61026203e-10   1.17210230e-10   8.53167852e-11
   6.21016942e-11   4.52035366e-11   3.29034456e-11   2.39502662e-11
   1.74332882e-11   1.26896100e-11   9.23670857e-12   6.72335754e-12
   4.89390092e-12   3.56224789e-12   2.59294380e-12   1.88739182e-12
   1.37382380e-12   1.00000000e-12]
[  1.00000000e-12   8.53167852e-13   7.27895384e-13   6.21016942e-13
   5.29831691e-13   4.52035366e-13   3.85662042e-13   3.29034456e-13
   2.80721620e-13   2.39502662e-13   2.04335972e-13   1.74332882e-13
   1.48735211e-13   1.26896100e-13   1.08263673e-13   9.23670857e-14
   7.88046282e-14   6.72335754e-14   5.73615251e-14   4.89390092e-14
   4.17531894e-14   3.56224789e-14   3.03919538e-14   2.59294380e-14
   2.21221629e-14   1.88739182e-14   1.61026203e-14   1.37382380e-14
   1.17210230e-14   1.00000000e-14]

In [67]:

    

xaxis3 = vary_x_samples3
yaxis3 = LogLikehood_wNoise_1e12(vary_x_samples3, sigma3)


fig3 = plt.figure()
plt.plot(xaxis3, yaxis3)
plt.xscale("log")
fig3.suptitle('Real-space LF vs C_3 values')
plt.xlabel('C_3 values, num=30, logspace from 1e-8 to 1e-12')
plt.ylabel('-2 ln likelihood L(T|C_3)')
plt.axvline(x =1.25873313e-09, color = 'g')  # a30 = 1.25873313e-09
# plt.plt.axvline(x = 4.79243265e-17 ) # a31 = 4.79243265e-17 
plt.axvline(x =5.86460794e-13, color = 'r')  # a32 = 5.86460794e-13
plt.axvline(x =7.27453196e-10, color = 'm')  # a33 = 7.27453196e-10
plt.axvline(x =4.5454678088e-09, color = 'y') # 7*\hat{C}_3 = 4.5454678088e-09
plt.axvline(x =6.49352544114e-10, color = 'k')  # \hat{C}_3 = 6.49352544114e-10
plt.grid()

print "vary_x_samples3 = np.logspace(-8, -12, num=30) #num = 30"
print "sigma3 = np.logspace(-12, -14, num=30)"
print "LogLikehood_wNoise_1e12(vary_x_samples3, sigma3)"
print "|a30|^2 is green, |a32|^2 is red, |a33|^2 is magenta, 7*\hat{C}_3 is yellow, \hat{C}_3 is black"


    

    




vary_x_samples3 = np.logspace(-8, -12, num=30) #num = 30
sigma3 = np.logspace(-12, -14, num=30)
LogLikehood_wNoise_1e12(vary_x_samples3, sigma3)
|a30|^2 is green, |a32|^2 is red, |a33|^2 is magenta, 7*\hat{C}_3 is yellow, \hat{C}_3 is black






    

In [68]:

    

print vary_x_samples3
print LogLikehood_wNoise_1e12(vary_x_samples3, sigma3)


    

    




[  1.00000000e-08   7.27895384e-09   5.29831691e-09   3.85662042e-09
   2.80721620e-09   2.04335972e-09   1.48735211e-09   1.08263673e-09
   7.88046282e-10   5.73615251e-10   4.17531894e-10   3.03919538e-10
   2.21221629e-10   1.61026203e-10   1.17210230e-10   8.53167852e-11
   6.21016942e-11   4.52035366e-11   3.29034456e-11   2.39502662e-11
   1.74332882e-11   1.26896100e-11   9.23670857e-12   6.72335754e-12
   4.89390092e-12   3.56224789e-12   2.59294380e-12   1.88739182e-12
   1.37382380e-12   1.00000000e-12]
[  7.13033005e+06   8.35364775e+06   9.78758542e+06   1.14683913e+07
   1.34385519e+07   1.57478658e+07   1.84547024e+07   2.16274760e+07
   2.53463756e+07   2.97053908e+07   3.48146872e+07   4.08033900e+07
   4.78228472e+07   5.60504544e+07   6.56941374e+07   7.69976068e+07
   9.02465162e+07   1.05775682e+08   1.23977542e+08   1.45312075e+08
   1.70318426e+08   1.99628525e+08   2.33983064e+08   2.74250209e+08
   3.21447543e+08   3.76767782e+08   4.41608915e+08   5.17609537e+08
   6.06690257e+08   7.11102239e+08]

In [69]:

    

#
# Plot noise sigma from e-20 to e-22
#
tempp = (1e6)*tempval # multiply CMB maps by 1e6

def LogLikehood_wNoise_1e20(param, sig):
    # param is our parameter, C_3
    Sij = param[:, None, None] * norm_matrix[1][None, :, :]
    newSij = (1e20)*Sij   # multiply S_ij by 1e12
    Nij = sig[:, None, None] * id_mat[None, :, :]
    newNij = (1e20)*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

#
# Plot noise sigma from e-20 to e-22
#
vary_x_samples12 = np.logspace(-8, -12, num=40) #num = 40
sigma12 = np.logspace(-20, -22, num=40)

xaxis12 = vary_x_samples12
yaxis12 = LogLikehood_wNoise_1e20(vary_x_samples12, sigma12)

fig12 = plt.figure()
plt.plot(xaxis12, yaxis12)
plt.xscale("log")
plt.xlabel('C_3 values, num=40, logspace from 1e-8 to 1e-12')
plt.ylabel('-2 ln likelihood L(T|C_3)')
plt.axvline(x =1.25873313e-09, color = 'g')  # a30 = 1.25873313e-09
# plt.plt.axvline(x = 4.79243265e-17 ) # a31 = 4.79243265e-17 
plt.axvline(x =5.86460794e-13, color = 'r')  # a32 = 5.86460794e-13
plt.axvline(x =7.27453196e-10, color = 'm')  # a33 = 7.27453196e-10
plt.axvline(x =4.5454678088e-09, color = 'y') # 7*\hat{C}_3 = 4.5454678088e-09
plt.axvline(x =6.49352544114e-10, color = 'k')  # \hat{C}_3 = 6.49352544114e-10
plt.grid()

print "vary_x_samples12 = np.logspace(-8, -12, num=40) #num = 40"
print "sigma12 = np.logspace(-20, -22, num=40)"
print "LogLikehood_wNoise_1e20(vary_x_samples12, sigma12)"
print "|a30|^2 is green, |a32|^2 is red, |a33|^2 is magenta, 7*\hat{C}_3 is yellow, \hat{C}_3 is black"


    

    




vary_x_samples12 = np.logspace(-8, -12, num=40) #num = 40
sigma12 = np.logspace(-20, -22, num=40)
LogLikehood_wNoise_1e20(vary_x_samples12, sigma12)
|a30|^2 is green, |a32|^2 is red, |a33|^2 is magenta, 7*\hat{C}_3 is yellow, \hat{C}_3 is black






    

In [70]:

    

print vary_x_samples12
print LogLikehood_wNoise_1e20(vary_x_samples12, sigma12)


    

    




[  1.00000000e-08   7.89652287e-09   6.23550734e-09   4.92388263e-09
   3.88815518e-09   3.07029063e-09   2.42446202e-09   1.91448198e-09
   1.51177507e-09   1.19377664e-09   9.42668455e-10   7.44380301e-10
   5.87801607e-10   4.64158883e-10   3.66524124e-10   2.89426612e-10
   2.28546386e-10   1.80472177e-10   1.42510267e-10   1.12533558e-10
   8.88623816e-11   7.01703829e-11   5.54102033e-11   4.37547938e-11
   3.45510729e-11   2.72833338e-11   2.15443469e-11   1.70125428e-11
   1.34339933e-11   1.06081836e-11   8.37677640e-12   6.61474064e-12
   5.22334507e-12   4.12462638e-12   3.25702066e-12   2.57191381e-12
   2.03091762e-12   1.60371874e-12   1.26638017e-12   1.00000000e-12]
[  7.27691114e+06   8.16795180e+06   9.16922170e+06   1.02974838e+07
   1.15663165e+07   1.29939553e+07   1.46026222e+07   1.64116138e+07
   1.84488568e+07   2.07409440e+07   2.33187451e+07   2.62201503e+07
   2.94881759e+07   3.31626926e+07   3.72983460e+07   4.19535599e+07
   4.71921318e+07   5.30870411e+07   5.97203796e+07   6.71863775e+07
   7.55881233e+07   8.50415290e+07   9.56802979e+07   1.07653711e+08
   1.21125729e+08   1.36287768e+08   1.53350921e+08   1.72551199e+08
   1.94158955e+08   2.18475972e+08   2.45840175e+08   2.76633026e+08
   3.11286032e+08   3.50281019e+08   3.94166775e+08   4.43550256e+08
   4.99122925e+08   5.61662545e+08   6.32040136e+08   7.11239239e+08]

In [71]:

    

#
# Plot noise sigma from e-21 to e-23
#
tempp = (1e6)*tempval # multiply CMB maps by 1e6

def LogLikehood_wNoise_1e21(param, sig):
    # param is our parameter, C_3
    Sij = param[:, None, None] * norm_matrix[1][None, :, :]
    newSij = (1e21)*Sij   # multiply S_ij by 1e12
    Nij = sig[:, None, None] * id_mat[None, :, :]
    newNij = (1e21)*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

#
# Plot noise sigma from e-21 to e-23
#
vary_x_samples125 = np.logspace(-8, -12, num=40) #num = 40
sigma125 = np.logspace(-21, -23, num=40)

xaxis125 = vary_x_samples125
yaxis125 = LogLikehood_wNoise_1e21(vary_x_samples125, sigma125)

fig125 = plt.figure()
plt.plot(xaxis125, yaxis125)
plt.xscale("log")
plt.xlabel('C_3 values, num=40, logspace from 1e-8 to 1e-12')
plt.ylabel('-2 ln likelihood L(T|C_3)')
plt.axvline(x =1.25873313e-09, color = 'g')  # a30 = 1.25873313e-09
# plt.plt.axvline(x = 4.79243265e-17 ) # a31 = 4.79243265e-17 
plt.axvline(x =5.86460794e-13, color = 'r')  # a32 = 5.86460794e-13
plt.axvline(x =7.27453196e-10, color = 'm')  # a33 = 7.27453196e-10
plt.axvline(x =4.5454678088e-09, color = 'y') # 7*\hat{C}_3 = 4.5454678088e-09
plt.axvline(x =6.49352544114e-10, color = 'k')  # \hat{C}_3 = 6.49352544114e-10
plt.grid()

print "vary_x_samples12 = np.logspace(-8, -12, num=40) #num = 40"
print "sigma12 = np.logspace(-21, -23, num=40)"
print "LogLikehood_wNoise_1e20(vary_x_samples12, sigma12)"
print "|a30|^2 is green, |a32|^2 is red, |a33|^2 is magenta, 7*\hat{C}_3 is yellow, \hat{C}_3 is black"


    

    




vary_x_samples12 = np.logspace(-8, -12, num=40) #num = 40
sigma12 = np.logspace(-21, -23, num=40)
LogLikehood_wNoise_1e20(vary_x_samples12, sigma12)
|a30|^2 is green, |a32|^2 is red, |a33|^2 is magenta, 7*\hat{C}_3 is yellow, \hat{C}_3 is black






    

In [72]:

    

print "vary_x_samples125 = np.logspace(-8, -12, num=40) #num = 40"
print "sigma125 = np.logspace(-21, -23, num=40)"
print "LogLikehood_wNoise_1e21(vary_x_samples125, sigma125)"
print vary_x_samples125
print LogLikehood_wNoise_1e21(vary_x_samples125, sigma125)


    

    




vary_x_samples125 = np.logspace(-8, -12, num=40) #num = 40
sigma125 = np.logspace(-21, -23, num=40)
LogLikehood_wNoise_1e21(vary_x_samples125, sigma125)
[  1.00000000e-08   7.89652287e-09   6.23550734e-09   4.92388263e-09
   3.88815518e-09   3.07029063e-09   2.42446202e-09   1.91448198e-09
   1.51177507e-09   1.19377664e-09   9.42668455e-10   7.44380301e-10
   5.87801607e-10   4.64158883e-10   3.66524124e-10   2.89426612e-10
   2.28546386e-10   1.80472177e-10   1.42510267e-10   1.12533558e-10
   8.88623816e-11   7.01703829e-11   5.54102033e-11   4.37547938e-11
   3.45510729e-11   2.72833338e-11   2.15443469e-11   1.70125428e-11
   1.34339933e-11   1.06081836e-11   8.37677640e-12   6.61474064e-12
   5.22334507e-12   4.12462638e-12   3.25702066e-12   2.57191381e-12
   2.03091762e-12   1.60371874e-12   1.26638017e-12   1.00000000e-12]
[  9.47758483e+06   1.01674782e+07   1.10456227e+07   1.20858191e+07
   1.32870283e+07   1.46310379e+07   1.61998771e+07   1.79770169e+07
   1.99725517e+07   2.22419458e+07   2.47805759e+07   2.76579653e+07
   3.09007610e+07   3.45488128e+07   3.86731641e+07   4.33261464e+07
   4.85426322e+07   5.44265135e+07   6.10528400e+07   6.85118369e+07
   7.69062033e+07   8.63395763e+07   9.70024716e+07   1.08958900e+08
   1.22435088e+08   1.37586869e+08   1.54644491e+08   1.73837852e+08
   1.95443578e+08   2.19756759e+08   2.47123994e+08   2.77918684e+08
   3.12551568e+08   3.51560133e+08   3.95433881e+08   4.44818141e+08
   5.00396903e+08   5.62937693e+08   6.33318209e+08   7.12514886e+08]

In [73]:

    

#
# Plot noise sigma from e-22 to e-24
#

def LogLikehood_wNoise_1e22(param, sig):
    # param is our parameter, C_3
    Sij = param[:, None, None] * norm_matrix[1][None, :, :]
    newSij = (1e22)*Sij   # multiply S_ij by 1e12
    Nij = sig[:, None, None] * id_mat[None, :, :]
    newNij = (1e22)*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_x_samples123 = np.logspace(-8, -12, num=40) #num = 40
sigma123 = np.logspace(-22, -24, num=40)

xaxis123 = vary_x_samples123
yaxis123 = LogLikehood_wNoise_1e22(vary_x_samples123, sigma123)

fig123 = plt.figure()
plt.plot(xaxis123, yaxis123)
plt.xscale("log")
plt.xlabel('C_3 values, num=40, logspace from 1e-8 to 1e-12')
plt.ylabel('-2 ln likelihood L(T|C_3)')
plt.axvline(x =1.25873313e-09, color = 'g')  # a30 = 1.25873313e-09
# plt.plt.axvline(x = 4.79243265e-17 ) # a31 = 4.79243265e-17 
plt.axvline(x =5.86460794e-13, color = 'r')  # a32 = 5.86460794e-13
plt.axvline(x =7.27453196e-10, color = 'm')  # a33 = 7.27453196e-10
plt.axvline(x =4.5454678088e-09, color = 'y') # 7*\hat{C}_3 = 4.5454678088e-09
plt.axvline(x =6.49352544114e-10, color = 'k')  # \hat{C}_3 = 6.49352544114e-10
plt.grid()

print "vary_x_samples123 = np.logspace(-8, -12, num=40) #num = 40"
print "sigma123 = np.logspace(-22, -24, num=40)"
print "LogLikehood_wNoise_1e22(vary_x_samples123, sigma123)"
print "|a30|^2 is green, |a32|^2 is red, |a33|^2 is magenta, 7*\hat{C}_3 is yellow, \hat{C}_3 is black"


    

    




vary_x_samples123 = np.logspace(-8, -12, num=40) #num = 40
sigma123 = np.logspace(-22, -24, num=40)
LogLikehood_wNoise_1e22(vary_x_samples123, sigma123)
|a30|^2 is green, |a32|^2 is red, |a33|^2 is magenta, 7*\hat{C}_3 is yellow, \hat{C}_3 is black






    

In [74]:

    

print vary_x_samples123
print LogLikehood_wNoise_1e22(vary_x_samples123, sigma123)


    

    




[  1.00000000e-08   7.89652287e-09   6.23550734e-09   4.92388263e-09
   3.88815518e-09   3.07029063e-09   2.42446202e-09   1.91448198e-09
   1.51177507e-09   1.19377664e-09   9.42668455e-10   7.44380301e-10
   5.87801607e-10   4.64158883e-10   3.66524124e-10   2.89426612e-10
   2.28546386e-10   1.80472177e-10   1.42510267e-10   1.12533558e-10
   8.88623816e-11   7.01703829e-11   5.54102033e-11   4.37547938e-11
   3.45510729e-11   2.72833338e-11   2.15443469e-11   1.70125428e-11
   1.34339933e-11   1.06081836e-11   8.37677640e-12   6.61474064e-12
   5.22334507e-12   4.12462638e-12   3.25702066e-12   2.57191381e-12
   2.03091762e-12   1.60371874e-12   1.26638017e-12   1.00000000e-12]
[  1.18423596e+06   5.05439716e+06   1.20679638e+06   3.80270233e+06
   2.09016375e+06   1.18209955e+06  -7.24125963e+05  -3.90564236e+06
  -9.91870220e+06  -2.29822343e+07  -5.63369055e+07  -3.44082279e+08
   2.13734532e+08   1.10217980e+08   9.28410346e+07   8.26604640e+07
   8.06174522e+07   8.18183755e+07   8.54764019e+07   8.98382791e+07
   9.79140866e+07   1.05974750e+08   1.15404844e+08   1.26577458e+08
   1.39150482e+08   1.53697442e+08   1.70526380e+08   1.89078287e+08
   2.10557047e+08   2.34471122e+08   2.61430158e+08   2.92212942e+08
   3.26649758e+08   3.65307711e+08   4.09121401e+08   4.58315512e+08
   5.13751370e+08   5.76402240e+08   6.46604825e+08   7.25788157e+08]

In [75]:

    

#
# Plot noise sigma from e-23 to e-25
#

def LogLikehood_wNoise_1e23(param, sig):
    # param is our parameter, C_3
    Sij = param[:, None, None] * norm_matrix[1][None, :, :]
    newSij = (1e23)*Sij   # multiply S_ij by 1e12
    Nij = sig[:, None, None] * id_mat[None, :, :]
    newNij = (1e23)*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_x_samples13 = np.logspace(-8, -12, num=40) #num = 40
sigma13 = np.logspace(-23, -25, num=40)

xaxis13 = vary_x_samples13
yaxis13 = LogLikehood_wNoise_1e23(vary_x_samples13, sigma13)

fig13 = plt.figure()
plt.plot(xaxis13, yaxis13)
plt.xscale("log")
plt.xlabel('C_3 values, num=40, logspace from 1e-8 to 1e-12')
plt.ylabel('-2 ln likelihood L(T|C_3)')
plt.axvline(x =1.25873313e-09, color = 'g')  # a30 = 1.25873313e-09
# plt.plt.axvline(x = 4.79243265e-17 ) # a31 = 4.79243265e-17 
plt.axvline(x =5.86460794e-13, color = 'r')  # a32 = 5.86460794e-13
plt.axvline(x =7.27453196e-10, color = 'm')  # a33 = 7.27453196e-10
plt.axvline(x =4.5454678088e-09, color = 'y') # 7*\hat{C}_3 = 4.5454678088e-09
plt.axvline(x =6.49352544114e-10, color = 'k')  # \hat{C}_3 = 6.49352544114e-10
plt.grid()

print "vary_x_samples13 = np.logspace(-8, -12, num=40) #num = 40"
print "sigma13 = np.logspace(-23, -25, num=40)"
print "LogLikehood_wNoise_1e23(vary_x_samples13, sigma13)"
print "|a30|^2 is green, |a32|^2 is red, |a33|^2 is magenta, 7*\hat{C}_3 is yellow, \hat{C}_3 is black"


    

    




vary_x_samples13 = np.logspace(-8, -12, num=40) #num = 40
sigma13 = np.logspace(-23, -25, num=40)
LogLikehood_wNoise_1e23(vary_x_samples13, sigma13)
|a30|^2 is green, |a32|^2 is red, |a33|^2 is magenta, 7*\hat{C}_3 is yellow, \hat{C}_3 is black






    

In [76]:

    

print "vary_x_samples13 = np.logspace(-8, -12, num=40) #num = 40"
print "sigma13 = np.logspace(-23, -25, num=40)"
print "LogLikehood_wNoise_1e23(vary_x_samples13, sigma13)"
print vary_x_samples13
print LogLikehood_wNoise_1e23(vary_x_samples13, sigma13)


    

    




vary_x_samples13 = np.logspace(-8, -12, num=40) #num = 40
sigma13 = np.logspace(-23, -25, num=40)
LogLikehood_wNoise_1e23(vary_x_samples13, sigma13)
[  1.00000000e-08   7.89652287e-09   6.23550734e-09   4.92388263e-09
   3.88815518e-09   3.07029063e-09   2.42446202e-09   1.91448198e-09
   1.51177507e-09   1.19377664e-09   9.42668455e-10   7.44380301e-10
   5.87801607e-10   4.64158883e-10   3.66524124e-10   2.89426612e-10
   2.28546386e-10   1.80472177e-10   1.42510267e-10   1.12533558e-10
   8.88623816e-11   7.01703829e-11   5.54102033e-11   4.37547938e-11
   3.45510729e-11   2.72833338e-11   2.15443469e-11   1.70125428e-11
   1.34339933e-11   1.06081836e-11   8.37677640e-12   6.61474064e-12
   5.22334507e-12   4.12462638e-12   3.25702066e-12   2.57191381e-12
   2.03091762e-12   1.60371874e-12   1.26638017e-12   1.00000000e-12]
[  3.77426185e+04  -4.71884682e+04   2.19525717e+06  -1.83104245e+06
   1.44077742e+06  -9.30100768e+05  -7.76671075e+05  -1.26936525e+05
   1.28253827e+07   9.68504741e+06  -8.76538959e+05  -1.47451062e+07
   5.16853388e+06   1.03687383e+07   2.63644475e+07  -5.73003279e+07
  -5.65150956e+06   4.77383307e+06  -1.02381080e+07   1.69372916e+07
  -2.72128212e+06   3.31160649e+07   3.90890261e+07   2.98301784e+07
   1.85441770e+07   3.69278683e+06  -2.01159266e+07  -6.85478208e+07
  -1.52581669e+08  -3.75382278e+08  -1.10801369e+09   6.65028850e+09
   1.37094621e+09   8.56374267e+08   8.52803804e+08   8.11743764e+08
   8.07777142e+08   8.37763508e+08   8.79143712e+08   9.46798929e+08]

In [77]:

    

#
# Plot noise sigma from e-24 to e-26
#

def LogLikehood_wNoise_1e24(param, sig):
    # param is our parameter, C_3
    Sij = param[:, None, None] * norm_matrix[1][None, :, :]
    newSij = (1e24)*Sij   # multiply S_ij by 1e12
    Nij = sig[:, None, None] * id_mat[None, :, :]
    newNij = (1e24)*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_x_samples14 = np.logspace(-8, -12, num=40) #num = 40
sigma14 = np.logspace(-24, -26, num=40)

xaxis14 = vary_x_samples14
yaxis14 = LogLikehood_wNoise_1e24(vary_x_samples14, sigma14)



fig14 = plt.figure()
plt.plot(xaxis14, yaxis14)
plt.xscale("log")
plt.xlabel('C_3 values, num=40, logspace from 1e-8 to 1e-12')
plt.ylabel('-2 ln likelihood L(T|C_3)')
plt.axvline(x =1.25873313e-09, color = 'g')  # a30 = 1.25873313e-09
# plt.plt.axvline(x = 4.79243265e-17 ) # a31 = 4.79243265e-17 
plt.axvline(x =5.86460794e-13, color = 'r')  # a32 = 5.86460794e-13
plt.axvline(x =7.27453196e-10, color = 'm')  # a33 = 7.27453196e-10
plt.axvline(x =4.5454678088e-09, color = 'y') # 7*\hat{C}_3 = 4.5454678088e-09
plt.axvline(x =6.49352544114e-10, color = 'k')  # \hat{C}_3 = 6.49352544114e-10
plt.grid()


print "vary_x_samples14 = np.logspace(-8, -12, num=40) #num = 40"
print "sigma14 = np.logspace(-24, -26, num=40)"
print "LogLikehood_wNoise(vary_x_samples14, sigma14)"
print "|a30|^2 is green, |a32|^2 is red, |a33|^2 is magenta, 7*\hat{C}_3 is yellow, \hat{C}_3 is black"


    

    




vary_x_samples14 = np.logspace(-8, -12, num=40) #num = 40
sigma14 = np.logspace(-24, -26, num=40)
LogLikehood_wNoise(vary_x_samples14, sigma14)
|a30|^2 is green, |a32|^2 is red, |a33|^2 is magenta, 7*\hat{C}_3 is yellow, \hat{C}_3 is black






    

In [78]:

    

print vary_x_samples14
print LogLikehood_wNoise_1e24(vary_x_samples14, sigma14)


    

    




[  1.00000000e-08   7.89652287e-09   6.23550734e-09   4.92388263e-09
   3.88815518e-09   3.07029063e-09   2.42446202e-09   1.91448198e-09
   1.51177507e-09   1.19377664e-09   9.42668455e-10   7.44380301e-10
   5.87801607e-10   4.64158883e-10   3.66524124e-10   2.89426612e-10
   2.28546386e-10   1.80472177e-10   1.42510267e-10   1.12533558e-10
   8.88623816e-11   7.01703829e-11   5.54102033e-11   4.37547938e-11
   3.45510729e-11   2.72833338e-11   2.15443469e-11   1.70125428e-11
   1.34339933e-11   1.06081836e-11   8.37677640e-12   6.61474064e-12
   5.22334507e-12   4.12462638e-12   3.25702066e-12   2.57191381e-12
   2.03091762e-12   1.60371874e-12   1.26638017e-12   1.00000000e-12]
[  4.14423224e+04  -2.68221082e+04  -4.24429210e+05  -2.10578450e+06
   2.47135612e+05  -6.06707744e+04   3.91217048e+05  -1.29760275e+06
   1.25524330e+05  -1.23160175e+05  -4.78897027e+05  -2.65982811e+06
   8.85507312e+05  -6.25246332e+04  -2.29375944e+06   2.61646465e+06
  -2.22803779e+06   6.65701789e+05  -2.13840424e+07  -2.77597765e+06
   2.69505537e+07   3.57189263e+06   2.94454370e+07  -9.01073786e+05
   1.92229256e+07  -7.72763940e+06  -9.09225749e+05  -3.74768615e+06
   6.73825231e+07  -5.52113110e+07   1.31605305e+08   3.62730958e+06
  -1.00644883e+06  -2.28428434e+07   3.15372407e+09  -1.91160005e+08
  -1.14152165e+08  -1.95995668e+08   1.54450820e+08   9.85681092e+08]

In [79]:

    

print hatC3  # anafast generated C_l for l = 3 from sky map
print 7*hatC3
print C3     # CAMB generated C_3, Boltzmann code output
print alm3_squared


    

    




6.49352544114e-10
4.5454678088e-09
5.88275e-10
[  1.25873313e-09   4.79243265e-17   5.86460794e-13   7.27453196e-10]

In [80]:

    

print tempval.shape


    

    




(3072,)

In [81]:

    

np.logspace(-8, -10, num=40)


    

    
Out[81]:





array([  1.00000000e-08,   8.88623816e-09,   7.89652287e-09,
         7.01703829e-09,   6.23550734e-09,   5.54102033e-09,
         4.92388263e-09,   4.37547938e-09,   3.88815518e-09,
         3.45510729e-09,   3.07029063e-09,   2.72833338e-09,
         2.42446202e-09,   2.15443469e-09,   1.91448198e-09,
         1.70125428e-09,   1.51177507e-09,   1.34339933e-09,
         1.19377664e-09,   1.06081836e-09,   9.42668455e-10,
         8.37677640e-10,   7.44380301e-10,   6.61474064e-10,
         5.87801607e-10,   5.22334507e-10,   4.64158883e-10,
         4.12462638e-10,   3.66524124e-10,   3.25702066e-10,
         2.89426612e-10,   2.57191381e-10,   2.28546386e-10,
         2.03091762e-10,   1.80472177e-10,   1.60371874e-10,
         1.42510267e-10,   1.26638017e-10,   1.12533558e-10,
         1.00000000e-10])

In [82]:

    

# Plot C_3 from e-08 to e-10
# Plot noise sigma from e-22 to e-23
#

def LogLikehood_wNoise_1e22(param, sig):
    # param is our parameter, C_3
    Sij = param[:, None, None] * norm_matrix[1][None, :, :]
    newSij = (1e22)*Sij   # multiply S_ij by 1e12
    Nij = sig[:, None, None] * id_mat[None, :, :]
    newNij = (1e22)*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_x_samples15 = np.logspace(-8, -10, num=40) #num = 40
sigma15 = np.logspace(-22, -23, num=40)

xaxis15 = vary_x_samples15
yaxis15 = LogLikehood_wNoise_1e22(vary_x_samples15, sigma15)

fig15 = plt.figure()
plt.plot(xaxis15, yaxis15)
plt.xscale("log")
plt.xlabel('C_3 values, num=40, logspace from 1e-8 to 1e-10')
plt.ylabel('-2 ln likelihood L(T|C_3)')
#plt.axvline(x =1.25873313e-09, color = 'g')  # a30 = 1.25873313e-09
# plt.plt.axvline(x = 4.79243265e-17 ) # a31 = 4.79243265e-17 
#plt.axvline(x =5.86460794e-13, color = 'r')  # a32 = 5.86460794e-13
#plt.axvline(x =7.27453196e-10, color = 'm')  # a33 = 7.27453196e-10
#plt.axvline(x =4.5454678088e-09, color = 'y') # 7*\hat{C}_3 = 4.5454678088e-09
#plt.axvline(x =6.49352544114e-10, color = 'k')  # \hat{C}_3 = 6.49352544114e-10
plt.grid()


print "|a30|^2 is green, |a32|^2 is red, |a33|^2 is magenta, 7*\hat{C}_3 is yellow, \hat{C}_3 is black"


    

    




|a30|^2 is green, |a32|^2 is red, |a33|^2 is magenta, 7*\hat{C}_3 is yellow, \hat{C}_3 is black






    

In [83]:

    

print vary_x_samples15
print LogLikehood_wNoise_1e22(vary_x_samples15, sigma15)


    

    




[  1.00000000e-08   8.88623816e-09   7.89652287e-09   7.01703829e-09
   6.23550734e-09   5.54102033e-09   4.92388263e-09   4.37547938e-09
   3.88815518e-09   3.45510729e-09   3.07029063e-09   2.72833338e-09
   2.42446202e-09   2.15443469e-09   1.91448198e-09   1.70125428e-09
   1.51177507e-09   1.34339933e-09   1.19377664e-09   1.06081836e-09
   9.42668455e-10   8.37677640e-10   7.44380301e-10   6.61474064e-10
   5.87801607e-10   5.22334507e-10   4.64158883e-10   4.12462638e-10
   3.66524124e-10   3.25702066e-10   2.89426612e-10   2.57191381e-10
   2.28546386e-10   2.03091762e-10   1.80472177e-10   1.60371874e-10
   1.42510267e-10   1.26638017e-10   1.12533558e-10   1.00000000e-10]
[  1.18423596e+06  -1.24653693e+04   5.05439716e+06   3.16566061e+06
   1.20679638e+06   8.44245200e+06   3.80270233e+06   3.26303568e+06
   2.09016375e+06   2.01331259e+06   1.18209955e+06   3.29437040e+05
  -7.24125963e+05  -1.92847638e+06  -3.90564236e+06  -6.64634227e+06
  -9.91870220e+06  -1.49515512e+07  -2.29822343e+07  -3.43676154e+07
  -5.63369055e+07  -1.16912692e+08  -3.44082279e+08   1.63685829e+08
   2.13734532e+08   1.38908195e+08   1.10217980e+08   8.53769658e+07
   9.28410346e+07   8.53345465e+07   8.26604640e+07   8.06413908e+07
   8.06174522e+07   8.06604458e+07   8.18183755e+07   8.33685071e+07
   8.54764019e+07   8.76337626e+07   8.98382791e+07   9.44712776e+07]

In [84]:

    

np.logspace(-8.3, -10, num=40)


    

    
Out[84]:





array([  5.01187234e-09,   4.53325612e-09,   4.10034607e-09,
         3.70877741e-09,   3.35460219e-09,   3.03424945e-09,
         2.74448928e-09,   2.48240018e-09,   2.24533968e-09,
         2.03091762e-09,   1.83697212e-09,   1.66154773e-09,
         1.50287575e-09,   1.35935639e-09,   1.22954263e-09,
         1.11212562e-09,   1.00592153e-09,   9.09859552e-10,
         8.22971159e-10,   7.44380301e-10,   6.73294595e-10,
         6.08997324e-10,   5.50840217e-10,   4.98236910e-10,
         4.50657034e-10,   4.07620869e-10,   3.68694506e-10,
         3.33485475e-10,   3.01638783e-10,   2.72833338e-10,
         2.46778711e-10,   2.23212211e-10,   2.01896228e-10,
         1.82615847e-10,   1.65176674e-10,   1.49402882e-10,
         1.35135431e-10,   1.22230472e-10,   1.10557891e-10,
         1.00000000e-10])

In [85]:

    

# Plot C_3 from 7.07945784e-09 to 1.00000000e-10
# Plot noise sigma from e-22 to e-23
#

def LogLikehood_wNoise_1e22(param, sig):
    # param is our parameter, C_3
    Sij = param[:, None, None] * norm_matrix[1][None, :, :]
    newSij = (1e22)*Sij   # multiply S_ij by 1e12
    Nij = sig[:, None, None] * id_mat[None, :, :]
    newNij = (1e22)*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_x_samples16 = np.logspace(-8.7, -10, num=40) #num = 40
sigma16 = np.logspace(-22, -24, num=40)

xaxis16 = vary_x_samples16
yaxis16 = LogLikehood_wNoise_1e22(vary_x_samples16, sigma16)

fig16 = plt.figure()
plt.plot(xaxis16, yaxis16)
plt.xscale("log")
plt.xlabel('C_3 values, num=40, logspace from 1e-9 to 1e-10')
plt.ylabel('-2 ln likelihood L(T|C_3)')
#plt.axvline(x =1.25873313e-09, color = 'g')  # a30 = 1.25873313e-09
# plt.plt.axvline(x = 4.79243265e-17 ) # a31 = 4.79243265e-17 
#plt.axvline(x =5.86460794e-13, color = 'r')  # a32 = 5.86460794e-13
#plt.axvline(x =7.27453196e-10, color = 'm')  # a33 = 7.27453196e-10
#plt.axvline(x =4.5454678088e-09, color = 'y') # 7*\hat{C}_3 = 4.5454678088e-09
#plt.axvline(x =6.49352544114e-10, color = 'k')  # \hat{C}_3 = 6.49352544114e-10
plt.grid()

vary_x_samples66 = np.logspace(-8.7, -11, num=40) #num = 40
sigma66 = np.logspace(-22, -23, num=40)

xaxis66 = vary_x_samples66
yaxis66 = LogLikehood_wNoise_1e22(vary_x_samples66, sigma66)

fig66 = plt.figure()
plt.plot(xaxis66, yaxis66)
plt.xscale("log")
plt.xlabel('C_3 values, num=40, logspace from 1e-9 to 1e-10')
plt.ylabel('-2 ln likelihood L(T|C_3)')
#plt.axvline(x =1.25873313e-09, color = 'g')  # a30 = 1.25873313e-09
# plt.plt.axvline(x = 4.79243265e-17 ) # a31 = 4.79243265e-17 
#plt.axvline(x =5.86460794e-13, color = 'r')  # a32 = 5.86460794e-13
#plt.axvline(x =7.27453196e-10, color = 'm')  # a33 = 7.27453196e-10
#plt.axvline(x =4.5454678088e-09, color = 'y') # 7*\hat{C}_3 = 4.5454678088e-09
#plt.axvline(x =6.49352544114e-10, color = 'k')  # \hat{C}_3 = 6.49352544114e-10
plt.grid()


print "|a30|^2 is green, |a32|^2 is red, |a33|^2 is magenta, 7*\hat{C}_3 is yellow, \hat{C}_3 is black"


    

    




|a30|^2 is green, |a32|^2 is red, |a33|^2 is magenta, 7*\hat{C}_3 is yellow, \hat{C}_3 is black






    













    

In [86]:

    

# Plot C_3 from 7.07945784e-09 to 1.00000000e-10
# Plot noise sigma from e-22 to e-23
#

def LogLikehood_wNoise_1e22(param, sig):
    # param is our parameter, C_3
    Sij = param[:, None, None] * norm_matrix[1][None, :, :]
    newSij = (1e22)*Sij   # multiply S_ij by 1e12
    Nij = sig[:, None, None] * id_mat[None, :, :]
    newNij = (1e22)*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_x_samples17 = np.logspace(-8.1, -10, num=40) #num = 40
sigma17 = np.logspace(-22, -23, num=40)

xaxis17 = vary_x_samples17
yaxis17 = LogLikehood_wNoise_1e22(vary_x_samples17, sigma17)

fig16 = plt.figure()
plt.plot(xaxis17, yaxis17)
plt.xscale("log")
plt.xlabel('C_3 values, num=40, logspace from 1e-9 to 1e-10')
plt.ylabel('-2 ln likelihood L(T|C_3)')
#plt.axvline(x =1.25873313e-09, color = 'g')  # a30 = 1.25873313e-09
# plt.plt.axvline(x = 4.79243265e-17 ) # a31 = 4.79243265e-17 
#plt.axvline(x =5.86460794e-13, color = 'r')  # a32 = 5.86460794e-13
#plt.axvline(x =7.27453196e-10, color = 'm')  # a33 = 7.27453196e-10
#plt.axvline(x =4.5454678088e-09, color = 'y') # 7*\hat{C}_3 = 4.5454678088e-09
#plt.axvline(x =6.49352544114e-10, color = 'k')  # \hat{C}_3 = 6.49352544114e-10
plt.grid()

vary_x_samples177 = np.logspace(-8.1, -10, num=40) #num = 40
sigma177 = np.logspace(-22, -24, num=40)

xaxis177 = vary_x_samples177
yaxis177 = LogLikehood_wNoise_1e22(vary_x_samples177, sigma177)

fig167 = plt.figure()
plt.plot(xaxis177, yaxis177)
plt.xscale("log")
plt.xlabel('C_3 values, num=40, logspace from 1e-9 to 1e-10')
plt.ylabel('-2 ln likelihood L(T|C_3)')
#plt.axvline(x =1.25873313e-09, color = 'g')  # a30 = 1.25873313e-09
# plt.plt.axvline(x = 4.79243265e-17 ) # a31 = 4.79243265e-17 
#plt.axvline(x =5.86460794e-13, color = 'r')  # a32 = 5.86460794e-13
#plt.axvline(x =7.27453196e-10, color = 'm')  # a33 = 7.27453196e-10
#plt.axvline(x =4.5454678088e-09, color = 'y') # 7*\hat{C}_3 = 4.5454678088e-09
#plt.axvline(x =6.49352544114e-10, color = 'k')  # \hat{C}_3 = 6.49352544114e-10
plt.grid()




print "|a30|^2 is green, |a32|^2 is red, |a33|^2 is magenta, 7*\hat{C}_3 is yellow, \hat{C}_3 is black"


    

    




|a30|^2 is green, |a32|^2 is red, |a33|^2 is magenta, 7*\hat{C}_3 is yellow, \hat{C}_3 is black






    













    

In [87]:

    

param = vary_x_samples17
sig = sigma17
Sij = param[:, None, None] * norm_matrix[1][None, :, :]
newSij = (1e22)*Sij   # multiply S_ij by 1e12
Nij = sig[:, None, None] * id_mat[None, :, :]
newNij = (1e22)*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]) ])
np.linalg.solve(Cij[i], tempp)
model_fit_terms + logdetC[1] + Npix2pi


    

    
Out[87]:





array([  3.38252302e+06,   2.57927165e+06,   4.09458156e+06,
         2.77756712e+06,   2.32039226e+06,   1.80336455e+06,
         1.38823728e+06,   9.72880913e+05,  -6.53682594e+03,
        -5.20254762e+05,  -1.68859911e+06,  -2.69025627e+06,
        -4.98145629e+06,  -7.55526647e+06,  -1.08993851e+07,
        -1.54105819e+07,  -2.21591599e+07,  -3.44419942e+07,
        -5.43932611e+07,  -9.46997897e+07,  -2.75873984e+08,
         3.51068588e+09,   2.26707975e+08,   1.43705216e+08,
         1.10924331e+08,   9.39441355e+07,   1.16056202e+08,
         8.40083422e+07,   7.96329150e+07,   7.82465573e+07,
         7.68246186e+07,   7.67974003e+07,   7.76128964e+07,
         7.85872039e+07,   8.00463879e+07,   8.22968098e+07,
         8.45469802e+07,   8.72974347e+07,   9.00461005e+07,
         9.44712776e+07])

In [88]:

    

print param.shape
print sig.shape
print Sij.shape
print Nij.shape
print norm_matrix[1].shape
print id_mat.shape
print model_fit_terms.shape
print np.linalg.solve(Cij[i], tempp).shape
print np.array([np.dot(tempp.T , np.linalg.solve(Cij[i], tempp) ) for i in range(Cij.shape[0]) ])
print np.array([np.dot(tempp.T , np.linalg.solve(Cij[i], tempp) ) for i in range(Cij.shape[0]) ]).shape


    

    




(40,)
(40,)
(40, 3072, 3072)
(40, 3072, 3072)
(3072, 3072)
(3072, 3072)
(40,)
(3072,)
[  3.36344611e+06   2.56030484e+06   4.07571941e+06   2.75886856e+06
   2.30182477e+06   1.78494492e+06   1.36997885e+06   9.54784509e+05
  -2.44858712e+04  -5.38035826e+05  -1.70621135e+06  -2.70769123e+06
  -4.99873407e+06  -7.57237650e+06  -1.09163194e+07  -1.54273424e+07
  -2.21757463e+07  -3.44584013e+07  -5.44094950e+07  -9.47158510e+07
  -2.75889849e+08   3.51067019e+09   2.26692457e+08   1.43689872e+08
   1.10909163e+08   9.39291477e+07   1.16041398e+08   8.39937136e+07
   7.96184653e+07   7.82322859e+07   7.68105259e+07   7.67834878e+07
   7.75991641e+07   7.85736545e+07   8.00330172e+07   8.22836224e+07
   8.45339736e+07   8.72846099e+07   9.00334591e+07   9.44588157e+07]
(40,)

In [89]:

    

invvv = np.array([np.linalg.solve(Cij[i], np.identity(3072)) for i in range(Cij.shape[0]) ])
#print invvv


    

In [90]:

    

# test of 
# model_fit_terms = np.array([np.dot(tempp.T , np.linalg.solve(Cij[i], tempp) ) for i in range(Cij.shape[0]) ])
# 
model_fit_terms22 = np.einsum('j,ijk,k->i',tempp, invvv, tempp)
print model_fit_terms22.shape
print model_fit_terms22


    

    




(40,)
[  3.36344611e+06   2.56030484e+06   4.07571941e+06   2.75886856e+06
   2.30182477e+06   1.78494492e+06   1.36997885e+06   9.54784509e+05
  -2.44858712e+04  -5.38035826e+05  -1.70621135e+06  -2.70769123e+06
  -4.99873407e+06  -7.57237650e+06  -1.09163194e+07  -1.54273424e+07
  -2.21757463e+07  -3.44584013e+07  -5.44094950e+07  -9.47158510e+07
  -2.75889849e+08   3.51067019e+09   2.26692457e+08   1.43689872e+08
   1.10909163e+08   9.39291477e+07   1.16041398e+08   8.39937136e+07
   7.96184653e+07   7.82322859e+07   7.68105259e+07   7.67834878e+07
   7.75991641e+07   7.85736545e+07   8.00330172e+07   8.22836224e+07
   8.45339736e+07   8.72846099e+07   9.00334591e+07   9.44588157e+07]

In [91]:

    

vary_x_samples18 = np.logspace(-8.1, -10, num=5) #num = 5
sigma18 = np.logspace(-22, -23, num=5)


    

In [92]:

    

def LogLikehood_wNoise_1e22_contour(param, sig):
    # param is our parameter, C_3
    Sij = param * norm_matrix[1][None, :, :]
    newSij = (1e22)*Sij   # multiply S_ij by 1e12
    Nij = sig * id_mat[None, :, :]
    newNij = (1e22)*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 [93]:

    

arr_plot1 = np.array([[LogLikehood_wNoise_1e22_contour(np.asarray(i), np.asarray(j)) for i in vary_x_samples18] for j in sigma18])


    

In [94]:

    

print arr_plot1


    

    




[[[  3.38252302e+06]
  [ -8.83962484e+08]
  [  9.04835818e+06]
  [  7.60737015e+06]
  [  7.27753934e+06]]

 [[  2.73999798e+06]
  [ -1.32768192e+06]
  [  2.28449406e+07]
  [  1.43468729e+07]
  [  1.31427689e+07]]

 [[ -1.17075118e+07]
  [  9.19545792e+06]
  [ -1.52722226e+08]
  [  2.92259872e+07]
  [  2.41025848e+07]]

 [[  2.76075737e+05]
  [  3.52431496e+07]
  [ -6.88962474e+05]
  [  7.77687113e+07]
  [  4.56945895e+07]]

 [[  2.06110002e+07]
  [  3.08522548e+06]
  [  1.39529430e+08]
  [ -2.36691023e+08]
  [  9.44712776e+07]]]

In [95]:

    

print arr_plot1.shape


    

    




(5, 5, 1)

In [96]:

    

print arr_plot1.reshape(5,5)


    

    




[[  3.38252302e+06  -8.83962484e+08   9.04835818e+06   7.60737015e+06
    7.27753934e+06]
 [  2.73999798e+06  -1.32768192e+06   2.28449406e+07   1.43468729e+07
    1.31427689e+07]
 [ -1.17075118e+07   9.19545792e+06  -1.52722226e+08   2.92259872e+07
    2.41025848e+07]
 [  2.76075737e+05   3.52431496e+07  -6.88962474e+05   7.77687113e+07
    4.56945895e+07]
 [  2.06110002e+07   3.08522548e+06   1.39529430e+08  -2.36691023e+08
    9.44712776e+07]]

In [97]:

    

vary_x_samples19 = np.logspace(-8.1, -10, num=10) #num = 10
sigma19 = np.logspace(-22, -23, num=10)

print vary_x_samples19
print sigma19


    

    




[  7.94328235e-09   4.88527357e-09   3.00453853e-09   1.84784980e-09
   1.13646367e-09   6.98947321e-10   4.29866235e-10   2.64376119e-10
   1.62596469e-10   1.00000000e-10]
[  1.00000000e-22   7.74263683e-23   5.99484250e-23   4.64158883e-23
   3.59381366e-23   2.78255940e-23   2.15443469e-23   1.66810054e-23
   1.29154967e-23   1.00000000e-23]

In [98]:

    

arr_plot2 = np.array([[LogLikehood_wNoise_1e22_contour(np.asarray(i), np.asarray(j)) for i in vary_x_samples19] for j in sigma19])


    

In [99]:

    

print arr_plot2.shape
#print arr_plot2


    

    




(10, 10, 1)

In [100]:

    

import pylab as pb    
import numpy as np
import matplotlib.pyplot as plt



zz = arr_plot2.reshape(10,10)

plt.figure()
CS = plt.contour(vary_x_samples19, sigma19, zz)
plt.clabel(CS, inline=1, fontsize=10)

pb.show()


    

    




/opt/local/Library/Frameworks/Python.framework/Versions/2.7/lib/python2.7/site-packages/matplotlib/text.py:52: UnicodeWarning: Unicode equal comparison failed to convert both arguments to Unicode - interpreting them as being unequal
  if rotation in ('horizontal', None):
/opt/local/Library/Frameworks/Python.framework/Versions/2.7/lib/python2.7/site-packages/matplotlib/text.py:54: UnicodeWarning: Unicode equal comparison failed to convert both arguments to Unicode - interpreting them as being unequal
  elif rotation == 'vertical':






    

In [101]:

    

vary_x_samples19 = np.logspace(-8.1, -10, num=10) #num = 10
sigma19 = np.logspace(-22, -23, num=10)

xxx3 = vary_x_samples19[:5]
yyy3 = sigma19[:5]

zzz3 = np.array([[LogLikehood_wNoise_1e22_contour(np.asarray(i), np.asarray(j)) for i in xxx3] for j in yyy3])


    

In [102]:

    

zzz3_reshaped = zzz3.reshape(5,5)

plt.figure()
CS = plt.contour(xxx3, yyy3, zzz3_reshaped)
plt.clabel(CS, inline=1, fontsize=10)

pb.show()

print xxx3
print yyy3
print zzz3_reshaped


    

    













    




[  7.94328235e-09   4.88527357e-09   3.00453853e-09   1.84784980e-09
   1.13646367e-09]
[  1.00000000e-22   7.74263683e-23   5.99484250e-23   4.64158883e-23
   3.59381366e-23]
[[  3.38252302e+06  -4.13805810e+05  -2.24728383e+07   1.60970967e+07
    9.91716392e+06]
 [  7.15800869e+07   1.88880967e+06  -5.06772408e+06   5.72241980e+07
    1.52417620e+07]
 [ -1.13324064e+06   5.18584833e+06  -3.60423743e+05  -2.89294053e+07
    2.83087564e+07]
 [  1.08698755e+07   9.47022324e+06   3.55129567e+06  -7.55526647e+06
    1.22262665e+08]
 [  6.85240223e+05   1.10988491e+05   2.90013913e+07  -1.82676007e+05
   -4.20677057e+07]]

In [103]:

    

vary_x_samples19 = np.logspace(-8.1, -10, num=10) #num = 10
sigma19 = np.logspace(-22, -23, num=10)

xxx3 = vary_x_samples19[:3]
yyy3 = sigma19[:3]

zzz3 = np.array([[LogLikehood_wNoise_1e22_contour(np.asarray(i), np.asarray(j)) for i in xxx3] for j in yyy3])


zzz3_reshaped = zzz3.reshape(3,3)

plt.figure()
CS = plt.contour(xxx3, yyy3, zzz3_reshaped)
plt.clabel(CS, inline=1, fontsize=10)

pb.show()

print xxx3
print yyy3
print zzz3


    

    













    




[  7.94328235e-09   4.88527357e-09   3.00453853e-09]
[  1.00000000e-22   7.74263683e-23   5.99484250e-23]
[[[  3382523.02181872]
  [  -413805.81039619]
  [-22472838.27094319]]

 [[ 71580086.91456631]
  [  1888809.66666395]
  [ -5067724.08335905]]

 [[ -1133240.63958953]
  [  5185848.33082371]
  [  -360423.74296772]]]

In [ ]:

    




    

In [104]:

    

vary_x_samples33 = np.logspace(-7, -11, num=10) #num = 10
sigma33 = np.logspace(-22, -23, num=10)

xxx4 = vary_x_samples33[:5]
yyy4 = sigma33[:5]

zzz4 = np.array([[LogLikehood_wNoise_1e22_contour(np.asarray(i), np.asarray(j)) for i in xxx4] for j in yyy4])


    

In [105]:

    

zzz4_reshaped = zzz4.reshape(5,5)

plt.figure()
CS = plt.contour(xxx4, yyy4, zzz4_reshaped)
plt.clabel(CS, inline=1, fontsize=10)

pb.show()

print xxx4
print yyy4
print zzz4_reshaped


    

    













    




[  1.00000000e-07   3.59381366e-08   1.29154967e-08   4.64158883e-09
   1.66810054e-09]
[  1.00000000e-22   7.74263683e-23   5.99484250e-23   4.64158883e-23
   3.59381366e-23]
[[ -1.86050500e+06   6.06722516e+04   9.70770908e+05  -9.44423979e+05
    1.36284179e+07]
 [  6.42066927e+05   1.01814473e+05   1.29852559e+07   1.86084207e+06
    3.07094786e+07]
 [ -1.01300674e+06  -1.17483263e+06   6.41138247e+05   4.24428106e+06
   -9.80482709e+07]
 [  9.73735015e+05   4.99434514e+05   7.49608213e+06  -6.51262903e+06
   -1.36906438e+07]
 [ -9.26202098e+04   6.57285318e+04   3.29834179e+06   5.26474689e+07
   -2.74797381e+06]]

In [106]:

    

vary_x_samples44 = np.logspace(-7, -11, num=20) #num = 10
sigma44 = np.logspace(-22, -23, num=20)

xxx44 = vary_x_samples44[:4]
yyy44 = sigma44[:4]

zzz44 = np.array([[LogLikehood_wNoise_1e22_contour(np.asarray(i), np.asarray(j)) for i in xxx44] for j in yyy44])

zzz44_reshaped = zzz44.reshape(4,4)

plt.figure()
CS = plt.contour(xxx44, yyy44, zzz44_reshaped)
plt.clabel(CS, inline=1, fontsize=10)

pb.show()

print xxx44
print yyy44
print zzz44


    

    













    




[  1.00000000e-07   6.15848211e-08   3.79269019e-08   2.33572147e-08]
[  1.00000000e-22   8.85866790e-23   7.84759970e-23   6.95192796e-23]
[[[ -1860504.99566932]
  [  3290193.90946934]
  [  2064995.42007263]
  [  3254836.54560861]]

 [[  -773032.63068625]
  [   -58864.01399538]
  [ 35606164.23995649]
  [  -789831.36262707]]

 [[  -205000.90176436]
  [   128232.99086233]
  [  -641279.23843855]
  [  -366591.81939573]]

 [[   659563.87681533]
  [  2644270.95149969]
  [   119527.55300234]
  [ 11207488.10158314]]]

In [107]:

    

vary_x_samples55 = np.logspace(-7, -11, num=10) #num = 10
sigma55 = np.logspace(-22, -23, num=10)

xxx55 = vary_x_samples55
yyy55 = sigma55

zzz55 = np.array([[LogLikehood_wNoise_1e22_contour(np.asarray(i), np.asarray(j)) for i in xxx55] for j in yyy55])

zzz55_reshaped = zzz55.reshape(10,10)

plt.figure()
CS = plt.contour(xxx55, yyy55, zzz55_reshaped)
plt.clabel(CS, inline=1, fontsize=10)

pb.show()

print xxx55
print yyy55
print zzz55_reshaped


    

    













    




[  1.00000000e-07   3.59381366e-08   1.29154967e-08   4.64158883e-09
   1.66810054e-09   5.99484250e-10   2.15443469e-10   7.74263683e-11
   2.78255940e-11   1.00000000e-11]
[  1.00000000e-22   7.74263683e-23   5.99484250e-23   4.64158883e-23
   3.59381366e-23   2.78255940e-23   2.15443469e-23   1.66810054e-23
   1.29154967e-23   1.00000000e-23]
[[ -1.86050500e+06   6.06722516e+04   9.70770908e+05  -9.44423979e+05
    1.36284179e+07   8.22914284e+06   7.46318215e+06   7.24302457e+06
    7.16984946e+06   7.14455816e+06]
 [  6.42066927e+05   1.01814473e+05   1.29852559e+07   1.86084207e+06
    3.07094786e+07   1.12047835e+07   9.77700764e+06   9.39141908e+06
    9.26939855e+06   9.22652819e+06]
 [ -1.01300674e+06  -1.17483263e+06   6.41138247e+05   4.24428106e+06
   -9.80482709e+07   1.57569506e+07   1.28582833e+07   1.21976748e+07
    1.19900599e+07   1.19190061e+07]
 [  9.73735015e+05   4.99434514e+05   7.49608213e+06  -6.51262903e+06
   -1.36906438e+07   2.30678192e+07   1.70701361e+07   1.58793759e+07
    1.55229219e+07   1.54028501e+07]
 [ -9.26202098e+04   6.57285318e+04   3.29834179e+06   5.26474689e+07
   -2.74797381e+06   3.79004658e+07   2.28754247e+07   2.07277253e+07
    2.01154048e+07   1.99128747e+07]
 [  2.29844240e+05  -1.13473659e+06  -9.56428803e+06   1.97573999e+06
    4.91148416e+06   8.46901770e+07   3.12437186e+07   2.71501120e+07
    2.60948439e+07   2.57550880e+07]
 [ -5.41098560e+05   3.68001911e+06   1.17211489e+06  -5.56241084e+06
    2.04265118e+07  -2.75048694e+08   4.39614817e+07   3.57467138e+07
    3.39032598e+07   3.33296317e+07]
 [ -5.41098560e+05   2.07088551e+05   1.28693584e+06  -4.02370881e+07
   -5.16609616e+07  -3.89498773e+07   6.44568742e+07   4.74080662e+07
    4.41462898e+07   4.31625656e+07]
 [ -1.02655674e+05   4.79298661e+06  -7.09262212e+06  -4.02919784e+07
    1.04777367e+07  -8.02286875e+06   1.06061435e+08   6.36218588e+07
    5.76345146e+07   5.59403443e+07]
 [ -1.02655674e+05  -6.59944516e+05   1.42668492e+07   2.23530763e+07
    2.68829322e+07   1.31267027e+07   2.34141930e+08   8.67897347e+07
    7.55155599e+07   7.25853360e+07]]

In [ ]:

    




    

In [108]:

    

vary_x_samples20 = np.logspace(-7, -10, num=10) #num = 10
sigma20 = np.logspace(-22, -24, num=10)

print vary_x_samples20
print sigma20


    

    




[  1.00000000e-07   4.64158883e-08   2.15443469e-08   1.00000000e-08
   4.64158883e-09   2.15443469e-09   1.00000000e-09   4.64158883e-10
   2.15443469e-10   1.00000000e-10]
[  1.00000000e-22   5.99484250e-23   3.59381366e-23   2.15443469e-23
   1.29154967e-23   7.74263683e-24   4.64158883e-24   2.78255940e-24
   1.66810054e-24   1.00000000e-24]

In [109]:

    

arr_plot3 = np.array([[LogLikehood_wNoise_1e22_contour(np.asarray(i), np.asarray(j)) for i in vary_x_samples20] for j in sigma20])


    

In [110]:

    

reshape_arr_plot3 = arr_plot3.reshape(10,10)

plt.figure()
CS = plt.contour(vary_x_samples19, sigma19, reshape_arr_plot3)
plt.clabel(CS, inline=1, fontsize=10)

pb.show()


    

In [ ]:

    




    

In [111]:

    

vary_x_samples555 = np.logspace(-7, -11, num=15) #num = 10
sigma555 = np.logspace(-22, -23, num=15)

xxx555 = vary_x_samples555
yyy555 = sigma555

zzz555 = np.array([[LogLikehood_wNoise_1e22_contour(np.asarray(i), np.asarray(j)) for i in xxx555] for j in yyy555])

zzz555_reshaped = zzz555.reshape(15,15)

plt.figure()
CS = plt.contour(xxx555, yyy555, zzz555_reshaped)
plt.clabel(CS, inline=1, fontsize=10)

pb.show()

#print xxx555
#print yyy555
#print zzz555_reshaped


    

In [112]:

    

vary_x_samples555 = np.logspace(-7, -11, num=15) #num = 10
sigma555 = np.logspace(-22, -23, num=15)

xxx555 = vary_x_samples555[:7]
yyy555 = sigma555[:7]

zzz555 = np.array([[LogLikehood_wNoise_1e22_contour(np.asarray(i), np.asarray(j)) for i in xxx555] for j in yyy555])

zzz555_reshaped = zzz555.reshape(7,7)

plt.figure()
CS = plt.contour(xxx555, yyy555, zzz555_reshaped)
plt.clabel(CS, inline=1, fontsize=10)

pb.show()

print xxx555
print yyy555
print zzz555_reshaped


    

    













    




[  1.00000000e-07   5.17947468e-08   2.68269580e-08   1.38949549e-08
   7.19685673e-09   3.72759372e-09   1.93069773e-09]
[  1.00000000e-22   8.48342898e-23   7.19685673e-23   6.10540230e-23
   5.17947468e-23   4.39397056e-23   3.72759372e-23]
[[ -1.86050500e+06  -1.59413385e+05  -2.05847336e+06  -3.86557263e+05
   -7.21232483e+05  -5.00436792e+06   1.80170077e+07]
 [ -7.73032631e+05   9.89162320e+06   3.10691485e+06  -8.81935229e+05
    4.17227705e+06  -2.04727034e+06   3.74029166e+07]
 [  6.42066927e+05  -2.83046752e+05   2.33689791e+06   3.62525394e+06
    3.87043161e+06   1.46757120e+05  -1.86151806e+08]
 [ -1.01300674e+06   3.71917635e+05   5.53871716e+06   9.89058660e+05
    8.96044966e+06   7.68643626e+05  -2.30029128e+07]
 [  1.04815317e+05  -4.56583957e+06   1.55632842e+08   2.63533418e+07
    3.29485063e+06   5.70130151e+06  -9.72598381e+06]
 [  9.73735015e+05   2.24369285e+06  -2.29466013e+06  -7.39178755e+05
    1.46867338e+07  -9.55792025e+05  -3.86206793e+06]
 [ -9.26202098e+04  -9.80420069e+05  -1.73234123e+06  -7.19105401e+06
   -2.30269671e+07   8.76078607e+06   4.44865051e+05]]

In [113]:

    

vary_x_samples88 = np.logspace(-6, -10, num=15) #num = 10
sigma88 = np.logspace(-22, -23, num=15)

xxx88 = vary_x_samples88
yyy88 = sigma88

zzz88 = np.array([[LogLikehood_wNoise_1e22_contour(np.asarray(i), np.asarray(j)) for i in xxx88] for j in yyy88])

zzz88_reshaped = zzz88.reshape(15, 15)

plt.figure()
CS = plt.contour(xxx88, yyy88, zzz88_reshaped)
plt.clabel(CS, inline=1, fontsize=10)

pb.show()

#print xxx88
#print yyy88
#print zzz88_reshaped