I have been having issues finding a good kernel for the emus, espeically the DS one. Gonna try a few things here.
In [22]:
import matplotlib
#matplotlib.use('Agg')
from matplotlib import pyplot as plt
%matplotlib inline
import seaborn as sns
sns.set()
In [23]:
from pearce.emulator import NashvilleHot
from itertools import product
from GPy.kern import *
import numpy as np
from os import path
from sys import argv
import h5py
In [24]:
training_file = '/u/ki/swmclau2/des/ds_hsab/PearceDsHSABCosmo.hdf5'
assert path.isfile(training_file)
test_file = '/u/ki/swmclau2/des/ds_hsab_test/PearceDsHSABCosmoTest.hdf5'
assert path.isfile(test_file)
#save_fname = argv[3]
In [25]:
with h5py.File(training_file, 'r') as f:
HOD_params = len(f.attrs['hod_param_names'])
In [26]:
cosmo_idx, hod_idx = 2,-1
cosmo_kernels = [Linear(input_dim=7, ARD=True), RBF(input_dim=7, ARD=True), Linear(input_dim=7, ARD=True) + RBF(input_dim=7, ARD=True), Linear(input_dim=7, ARD=True) + Matern32(input_dim=7, ARD=True), \
Matern32(input_dim=7, ARD=True)+RBF(input_dim=7, ARD=True) + Bias(input_dim=7)]
HOD_kernels = [ Matern32(input_dim=HOD_params, ARD=True), RBF(input_dim=HOD_params, ARD=True) + Linear(input_dim=HOD_params, ARD=True), Matern32(input_dim=HOD_params, ARD=True)+RBF(input_dim=HOD_params, ARD=True) + Bias(input_dim=HOD_params)\
, RBF(input_dim=HOD_params, ARD=True) ]
#k = (cosmo_kernels[3], HOD_kernels[0])
k = (cosmo_kernels[cosmo_idx], HOD_kernels[hod_idx])
In [27]:
print k
In [28]:
fixed_params = {'z':0.0}
# best for wp: either RBF & Matern32 or Linear + RBF & RBF & Matern32
# best for DS: Linear + RBF, Matern32
#kernels = product(cosmo_kernels, HOD_kernels)
#for k in kernels:
hyperparams = {'kernel': k , \
'optimize': True}
In [53]:
emu = NashvilleHot(training_file, hyperparams=hyperparams,fixed_params = fixed_params, downsample_factor = 0.01)
In [56]:
pred_y, data_y = emu.goodness_of_fit(test_file, statistic = None)#, downsample_factor = 0.1)
In [57]:
plt.plot(emu.scale_bin_centers, ((10**pred_y - 10**data_y)/(10**data_y)).mean(axis=1), label= 'Bias')
plt.plot(emu.scale_bin_centers, (np.abs(10**pred_y - 10**data_y)/(10**data_y)).mean(axis =1), label = 'Acc')
plt.plot(emu.scale_bin_centers, np.zeros_like(emu.scale_bin_centers))
plt.xscale('log')
plt.legend(loc='best')
plt.ylim([-0.05, 0.2])
plt.show()
# average over realizations
In [32]:
pred_y_rs= pred_y.reshape((len(emu.scale_bin_centers),5,7, -1), order = 'F')[:,0,:,:]
data_y_rs= data_y.reshape((len(emu.scale_bin_centers),5,7, -1), order = 'F').mean(axis = 1)
R = (10**pred_y_rs - 10**data_y_rs).reshape((18,-1), order = 'F')
In [33]:
cov = R.dot(R.T)/(R.shape[1]-1)
print k
print 'Yerr', np.sqrt(np.diag(cov))/(10**data_y.mean(axis=1))
print '*'*10
#np.save(save_fname, cov)
In [34]:
fname = '/u/ki/swmclau2/des/PearceMCMC/HOD_wp_ds_rmin_None_HOD.hdf5'
f = h5py.File(fname, 'r')
In [35]:
#if obs == 'wp':
# true_data = f['data'][:len(r_bins)-1]
#else:
mcmc_data = f['data'][18:]
#if obs == 'wp':
# cov = f['cov'][:len(r_bins)-1][:, :len(r_bins)-1]
#else:
cov = f['cov'][18:][:, 18:]
yerr = np.sqrt(np.diag(cov))
In [36]:
sim_info = eval(f.attrs['sim'])
cosmo_params = {'simname': sim_info['simname'],
'boxno': sim_info['sim_hps']['boxno'],\
'realization': sim_info['sim_hps']['realization'],
'scale_factors':[sim_info['scale_factor']],\
'system': sim_info['sim_hps']['system']}
hod_params = sim_info['hod_params']
f.close()
In [37]:
sim_info
Out[37]:
In [38]:
from pearce.mocks import cat_dict
boxno = 3
realization = 2
cat = cat_dict['testbox'](boxno = boxno, realization = realization )#construct the specified catalog!
In [39]:
hod_params = {'alpha': 1.083,
'conc_gal_bias': 1.0,
'logM0': 13.2,
'logM1': 14.2,
'logMmin': 13.06094217689994,
'sigma_logM': 0.2}
In [40]:
hod_idx = 64
In [41]:
with h5py.File(test_file, 'r') as f:
hod_param_names = f.attrs['hod_param_names']
hod_param_vals = f.attrs['hod_param_vals'][hod_idx]
cosmo_param_names = f.attrs['cosmo_param_names']
cosmo_param_vals = f.attrs['cosmo_param_vals'][boxno*5+realization]
true_data = f['cosmo_no_%02d'%(boxno*5+realization)]['a_1.000']['obs'][hod_idx]
hod_params = dict(zip(hod_param_names, hod_param_vals))
#hod_params.update(dict(zip(cosmo_param_names, cosmo_param_vals)))
In [42]:
cpv = cat._get_cosmo_param_names_vals()
cat_val_dict = {key: val for key, val in zip(cpv[0], cpv[1])}
In [43]:
true_param_dict = cat_val_dict.copy()
for hp, hv in hod_params.iteritems():
if hp == 'logMmin':
continue
true_param_dict[hp] = hv
true_pred = emu.emulate_wrt_r(true_param_dict)[0]
In [44]:
rbc = emu.scale_bin_centers
In [45]:
N = 10
cmap = sns.color_palette("BrBG_d", N)
In [46]:
fig = plt.figure(figsize=(10,7))
plt.plot(rbc, 10**true_pred/10**true_data,label = 'Emu at Truth', color ='k')
#plt.plot(rbc, 10**pop_xi.mean(axis = 0), label = 'Sim' )
#plt.errorbar(rbc, np.ones_like(true_data), yerr=yerr/true_data, label = 'Data')
plt.plot(rbc, np.ones_like(true_data), label = 'Data')
plt.xscale('log')
plt.legend(loc='best')
plt.show();
In [47]:
fig = plt.figure(figsize=(10,7))
varied_pname = 'omch2'
lower, upper = emu.get_param_bounds(varied_pname)
i = 0
for c, val in zip(cmap, np.linspace(lower, upper, N) ):
print i, val
i+=1
param_dict = true_param_dict.copy()
param_dict[varied_pname] = val
pred = emu.emulate_wrt_r(param_dict)[0]
plt.plot(rbc, (10**pred-10**true_data[-len(emu.scale_bin_centers):])/10**true_data[-len(emu.scale_bin_centers):],\
alpha = 0.5,label = val, color =c)
pred = emu.emulate_wrt_r(true_param_dict)[0]
plt.plot(rbc, (10**pred-10**true_data[-len(emu.scale_bin_centers):])/10**true_data[-len(emu.scale_bin_centers):], label = 'Truth', color = 'k')
#plt.errorbar(rbc, np.zeros_like(true_data[-len(emu.scale_bin_centers):]), yerr=yerr/true_data[-len(emu.scale_bin_centers):], label = 'Data')
#plt.loglog()
plt.xscale('log')
plt.legend(loc='best')
plt.show();
In [ ]:
In [ ]:
In [ ]: