# ---
# jupyter:
#   jupytext:
#     text_representation:
#       extension: .py
#       format_name: light
#       format_version: '1.5'
#       jupytext_version: 1.14.4
#   kernelspec:
#     display_name: Python 2
#     language: python
#     name: python2
# ---

# <h1>Metodo de Maxima Verossimilhança usando Espectro de Potência</h1>
# <h3>Aluno:<NAME>, Astroinformática 2021<h3>

# Nessa etapa utilizaremos uma simulação de um catalogo de objetos gerados a partir de um modelo de universo baseado nos parametros obtidos pelos resultados de 2018 do Planck. 

# importando o modulo classy 
# %matplotlib inline
import matplotlib.pyplot as plt
import numpy as np
from classy import Class 
from math import pi
from nbodykit.lab import *
from nbodykit import style, setup_logging
from nbodykit.lab import cosmology
# %config InlineBackend.figure_format = 'retina'
from scipy.interpolate import InterpolatedUnivariateSpline
plt.style.use(style.notebook)
setup_logging()

# +
redshift = 0.01
cosmo = cosmology.Cosmology(h=0.67556, T0_cmb=2.7255, Omega0_b = 0.0483,Omega0_cdm = 0.2638,n_s = 0.9619)
Plin = cosmology.LinearPower(cosmo, redshift, transfer='EisensteinHu')
b1 = 2.0

cat = LogNormalCatalog(Plin=Plin, nbar=3e-4, BoxSize=1380., Nmesh=256, bias=b1, seed=42)
# -

# add RSD
line_of_sight = [0,0,1]
cat['RSDPosition'] = cat['Position'] + cat['VelocityOffset'] * line_of_sight

# convert to a MeshSource, using TSC interpolation on 256^3 mesh
mesh = cat.to_mesh(window='tsc', Nmesh=256, compensated=True, position='RSDPosition')

plt.imshow(mesh.preview(axes=[0,1], Nmesh=200))

r = FFTPower(mesh, mode='1d', dk=0.005, kmin=0.01)

Pk = r.power

# +
# print the shot noise subtracted P(k)
k = np.logspace(-2, 0, 100)
plt.loglog(Pk['k'], Pk['power'].real - Pk.attrs['shotnoise'])
plt.loglog(k, Plin(k), c='k')

# format the axes
plt.xlabel(r"$k$ [$h \ \mathrm{Mpc}^{-1}$]")
plt.ylabel(r"$P(k)$ [$h^{-3}\mathrm{Mpc}^3$]")
plt.xlim(0.01, 0.6)
# -

power = Pk['power'].real - Pk.attrs['shotnoise']
Knovo = Pk['k']
print(type(power))
print(type(Knovo))
import scipy 

Pknovo = scipy.interpolate.interp1d(Knovo, power)

# +
# print the shot noise subtracted P(k)
k = np.logspace(-2, 0, 100)
plt.loglog(Pk['k'], Pk['power'].real - Pk.attrs['shotnoise'] ,c='k')
plt.loglog(k, Plin(k), c='b')
#plt.loglog(Knovo, Pknovo(Knovo), c='r')

# format the axes
plt.xlabel(r"$k$ [$h \ \mathrm{Mpc}^{-1}$]")
plt.ylabel(r"$P(k)$ [$h^{-3}\mathrm{Mpc}^3$]")
plt.xlim(0.01, 0.6)


# -

def NewSigma8(Pk, k,r=8, kmin=1e-5, kmax=1e1):
    import mcfit
    from scipy.interpolate import InterpolatedUnivariateSpline as spline

    k = numpy.logspace(numpy.log10(kmin), numpy.log10(kmax), len(Pk))
    #Pk = self(k)
    R, sigmasq = mcfit.TophatVar(k, lowring=True)(Pk, extrap=True)
    return spline(R, sigmasq)(r)**0.5


Nbins = 20 
SigmaMesh=NewSigma8(power,Knovo)
print(SigmaMesh)
len(power)

OmegaB = np.linspace(0.011, 0.06, Nbins)
logL = np.empty(Nbins)
sigmaSimulados = np.empty(Nbins)
for m in range(Nbins):
    cosmo = cosmology.Cosmology(h=0.67556, T0_cmb=2.7255, Omega0_b = OmegaB[m],Omega0_cdm = 0.2638,n_s = 0.9619)
    Plin = cosmology.LinearPower(cosmo, redshift, transfer='CLASS')
    sigmaSimulados[m] = Plin.sigma_r(r=8)
    ErroSigma = Plin.sigma_r(r=8)/10
    logL[m] =  -np.sum((0.5 * ((Plin.sigma_r(r=8) - SigmaMesh ) / ErroSigma) ** 2)) 

#Buscando qual foi o valor de OmegaB que determinou o máximo da Likelihood 
loc = np.where(logL == np.max(logL))
#print(loc)
#print(sigmaSimulados[loc])
print("O valor de Omega Barions que melhor ajusta os dados é {}".format(OmegaB[loc]))
print("E o besfit de sigma8 é {}".format(sigmaSimulados[loc]))
plt.axvline(x = sigmaSimulados[loc], color = 'b', label = 'Besfit Sigma8 ')
plt.plot(sigmaSimulados, logL, c='k')
plt.xlabel(r'$\sigma_{8}$')
plt.ylabel(r'$\mathcal{L}$')
plt.legend()

#Usando a mesma ídeia para restringir H0 utilizando um Omega0_b maior (0.06)
hSample = np.linspace(0.4, 0.9, Nbins)
logL = np.empty(Nbins)
sigmaSimulados = np.empty(Nbins)
for m in range(Nbins):
    cosmo = cosmology.Cosmology(h=hSample[m], T0_cmb=2.7255, Omega0_b = 0.0483,Omega0_cdm = 0.2638,n_s = 0.9619)
    Plin = cosmology.LinearPower(cosmo, redshift, transfer='CLASS')
    sigmaSimulados[m] = Plin.sigma_r(r=8)
    ErroSigma = Plin.sigma_r(r=8)/10
    logL[m] = - np.sum(0.5 * ((Plin.sigma_r(r=8) - SigmaMesh ) / ErroSigma) ** 2)  

#Buscando qual foi o valor de h que determinou o máximo da Likelihood 
loc = np.where(logL == np.max(logL))
#print(loc)
#print(sigmaSimulados[loc])
print("O valor de h que melhor ajusta os dados é {}".format(hSample[loc]))
print("E o besfit de sigma8 é {}".format(sigmaSimulados[loc]))
plt.axvline(x = sigmaSimulados[loc], color = 'b', label = 'Besfit Sigma8 ')
plt.plot(sigmaSimulados, logL, c='k')
plt.xlabel(r'$\sigma_{8}$')
plt.ylabel(r'$\mathcal{L}$')
plt.legend()

# +
#Tentando fitar Omega_b para o mesmo sigma8 obtido através da simulação porém utilizando uma cosmologia com H0 = 73
# -

OmegaB = np.linspace(0.011, 0.06, Nbins)
logL = np.empty(Nbins)
sigmaSimulados = np.empty(Nbins)
for m in range(Nbins):
    cosmo = cosmology.Cosmology(h=0.73, T0_cmb=2.7255, Omega0_b = OmegaB[m],Omega0_cdm = 0.2638,n_s = 0.9619)
    Plin = cosmology.LinearPower(cosmo, redshift, transfer='CLASS')
    sigmaSimulados[m] = Plin.sigma_r(r=8)
    ErroSigma = Plin.sigma_r(r=8)/10
    logL[m] =  -np.sum((0.5 * ((Plin.sigma_r(r=8) - SigmaMesh ) / ErroSigma) ** 2)) 

#Buscando qual foi o valor de OmegaB que determinou o máximo da Likelihood 
loc = np.where(logL == np.max(logL))
#print(loc)
#print(sigmaSimulados[loc])
print("O valor de Omega Barions que melhor ajusta os dados é {}".format(OmegaB[loc]))
print("E o besfit de sigma8 é {}".format(sigmaSimulados[loc]))
plt.axvline(x = sigmaSimulados[loc], color = 'b', label = 'Besfit Sigma8 ')
plt.plot(sigmaSimulados, logL, c='k')
plt.xlabel(r'$\sigma_{8}$')
plt.ylabel(r'$\mathcal{L}$')
plt.legend()



StarcoderdataJupyterScriptsDedupFiltered