Cosmo MCMC is landing on biased HOD + Cosmology. I'm gonna first look at the code to make preds, and compare to the emulator.
In [1]:
from pearce.emulator import OriginalRecipe, ExtraCrispy, SpicyBuffalo
from pearce.mocks import cat_dict
import numpy as np
from os import path
In [2]:
import matplotlib
#matplotlib.use('Agg')
from matplotlib import pyplot as plt
%matplotlib inline
import seaborn as sns
sns.set()
In [3]:
#training_file = '/scratch/users/swmclau2/xi_zheng07_cosmo_lowmsat/PearceRedMagicXiCosmoFixedNd.hdf5'
training_file = '/u/ki/swmclau2/des/PearceRedMagicXiCosmoFixedNdLowMsat.hdf5'
#test_file = '/scratch/users/swmclau2/xi_zheng07_cosmo_test_lowmsat/PearceRedMagicXiCosmoFixedNd_Test.hdf5'
test_file = '/u/ki/swmclau2/des/PearceRedMagicXiCosmoFixedNdLowMsatTest.hdf5'
#test_file = '/u/ki/swmclau2/des/xi_cosmo_tester/PearceRedMagicXiCosmoFixedNd_test.hdf5'
em_method = 'gp'
split_method = 'random'
In [4]:
a = 1.0
z = 1.0/a - 1.0
In [5]:
fixed_params = {'z':z}#, 'cosmo': 3}#, 'r':0.53882047}
In [6]:
np.random.seed(0)
emu = SpicyBuffalo(training_file, method = em_method, fixed_params=fixed_params,
custom_mean_function = 'linear', downsample_factor = 0.1)
In [7]:
emu.downsample_x[0]
Out[7]:
In [8]:
emu.get_param_names()
Out[8]:
In [9]:
v = np.ones_like(emu._emulators[0].get_parameter_vector())*12.0
In [10]:
if hasattr(emu, "_emulator"):
emu._emulator.set_parameter_vector(v)
emu._emulator.recompute()
else:
for _emulator in emu._emulators:
_emulator.set_parameter_vector(v)
_emulator.recompute()
In [11]:
emu._emulators[0].get_parameter_vector()
Out[11]:
In [12]:
gof = emu.goodness_of_fit(test_file, statistic='log_frac')
In [13]:
for g in gof:
print np.mean(g), np.median(g)
In [14]:
params = {}
for pname in emu.get_param_names():
if pname == 'r':
continue
low, high = emu.get_param_bounds(pname)
params[pname] = np.random.uniform(low, high)
print params
In [15]:
fixed_params = {}#'f_c':1.0}#,'logM1': 13.8 }# 'z':0.0}
cosmo_params = {'simname':'testbox', 'boxno': 3, 'realization': 0, 'scale_factors':[1.0], 'system': 'ki-ls'}
cat = cat_dict[cosmo_params['simname']](**cosmo_params)#construct the specified catalog!
In [16]:
cat.load(1.0, HOD='zheng07')
In [17]:
print len(cat.halocat.halo_table)
In [18]:
cat._get_cosmo_param_names_vals()[1]
Out[18]:
In [19]:
test_point_idx = 1
test_point_dict = dict(zip(emu.get_param_names(), emu.x[0][test_point_idx]*emu._x_std[0]+emu._x_mean[0]))
In [20]:
#hod_param_names = ['logM0', 'sigma_logM', 'logM1', 'alpha']
emulation_point = [('logM0', 14.0), ('sigma_logM', 0.2),
('alpha', 1.083),('logM1', 13.7)]#, ('logMmin', 12.233)]
#em_params = {key:test_point_dict[key] for key in hod_param_names}
#em_params = dict(zip(hod_param_names, x_point))
em_params = dict(emulation_point)
em_params.update(fixed_params)
In [21]:
from scipy.optimize import minimize_scalar
def add_logMmin(hod_params, cat, nd = 1e-4):
"""
In the fixed number density case, find the logMmin value that will match the nd given hod_params
:param: hod_params:
The other parameters besides logMmin
:param cat:
the catalog in question
:return:
None. hod_params will have logMmin added to it.
"""
hod_params['logMmin'] = 13.0 #initial guess
#cat.populate(hod_params) #may be overkill, but will ensure params are written everywhere
def func(logMmin, hod_params):
hod_params.update({'logMmin':logMmin})
calc_nd = cat.calc_analytic_nd(hod_params)
#print logMmin, calc_nd
return (calc_nd - nd)**2
res = minimize_scalar(func, bounds = (12, 16), args = (hod_params,), options = {'maxiter':100}, method = 'Bounded' )
# assuming this doens't fail
hod_params['logMmin'] = res.x
In [22]:
add_logMmin(em_params, cat)
In [23]:
r_bins = np.logspace(-1.1, 1.6, 19)
rpoints = (r_bins[1:]+r_bins[:-1])/2.0
In [24]:
em_params['logMmin']
Out[24]:
In [25]:
print em_params
In [26]:
# get cosmo params
try:
del em_params['logMmin']
except KeyError:
pass
cpv = cat._get_cosmo_param_names_vals()
cosmo_param_dict = {key: val for key, val in zip(cpv[0], cpv[1])}
em_params.update( cosmo_param_dict)
In [27]:
print em_params
In [28]:
y_pred = emu.emulate_wrt_r(em_params)[0]
In [29]:
y_calc = np.loadtxt('/u/ki/swmclau2/Git/pearce/bin/mcmc/xi_gg_true.npy')
In [30]:
fig = plt.figure(figsize = (12, 5))
plt.subplot(121)
#plt.plot(rpoints, y_calc, label = 'Sim')
plt.plot(rpoints, y_calc, label = 'Sim')
plt.plot(rpoints, y_pred, label = 'Emu')
plt.xscale('log')
plt.legend(loc = 'best')
plt.ylabel(r'$\xi_{gg}(r)$')
plt.xlabel(r'$r$ [Mpc]')
plt.subplot(122)
#plt.plot(rpoints, y_calc/y_point, label = 'Sim')
plt.plot(rpoints, (10**y_pred)/(10**y_calc), label = 'Emu/Sim')
plt.legend(loc = 'best')
plt.xlabel(r'$r$ [Mpc]')
plt.xscale('log')
plt.show()
In [112]:
print y_pred/y_calc
In [69]:
fixed_params = {'z':z, 'cosmo': 3}#, 'r':0.53882047}
train_x, train_y, _, info = emu.get_data(test_file, fixed_params, None)#, skip_nans = False)
In [70]:
info
Out[70]:
In [71]:
cpv = cat._get_cosmo_param_names_vals()
cosmo_params = dict(zip(cpv[0], cpv[1]))
In [72]:
for idx in xrange(100):
x_point = train_x[idx*emu.n_bins, :-1]
y_point = train_y[idx*emu.n_bins:(idx+1)*emu.n_bins]
pop_params = dict(zip(info['ordered_params'].keys(), x_point))
#add_logMmin(pop_params, cat)
print pop_params
#_xi_vals = []
#for i in xrange(10):
# cat.populate(pop_params, min_ptcl=100)
# _xi_vals.append(cat.calc_xi(r_bins))
#xi_vals = np.log10(np.array(_xi_vals))
#y_calc = xi_vals.mean(axis = 0)
pop_params.update(cosmo_params)
#del pop_params['logMmin']
y_pred = emu.emulate_wrt_r(pop_params)[0]
fig = plt.figure(figsize = (14, 6))
plt.subplot(121)
#plt.plot(rpoints, y_calc, label = 'Sim')
plt.plot(rpoints, y_point, label = 'Data')
plt.plot(rpoints, y_pred, label = 'Emu')
plt.xscale('log')
plt.legend(loc = 'best')
plt.subplot(122)
#plt.plot(rpoints, y_calc/y_point, label = 'Sim')
plt.plot(rpoints, y_pred/y_point, label = 'Emu/True')
plt.legend(loc = 'best')
plt.xscale('log')
plt.show()
In [ ]:
In [ ]:
In [ ]: