In [4]:
import pylab
import numpy as np
%pylab inline
from halotools.empirical_models import Zheng07 , model_defaults
from astropy.table import Table
import plotting
from plotting import prettyplot
prettyplot()
In [587]:
model = Zheng07(threshold = -21)
data = np.loadtxt("richness.dat")
model.param_dict
Out[587]:
In [5]:
def richness(group_id):
gals = Table()
gals['groupid'] = group_id
gals['dummy'] = 1
grouped_table = gals.group_by('groupid')
grp_richness = grouped_table['dummy'].groups.aggregate(np.sum)
return grp_richness
In [232]:
model2 = Zheng07(threshold=-21)
model2.param_dict["logM0"] = 12.1
In [233]:
model.populate_mock()
model2.populate_mock()
In [234]:
g2 = model2.mock.compute_fof_group_ids()
g1 = model.mock.compute_fof_group_ids()
In [589]:
bins = np.arange(1,61,1)
In [590]:
np.histogram(data,bins)
Out[590]:
In [221]:
x = plt.hist(data , bins , normed="True")[0]
In [269]:
y = np.histogram(richness(g2),bins , normed="True")[0]
yy = plt.hist(richness(g2),bins , normed="True")[0]
print yy.shape
z = plt.hist(richness(g1),bins , normed="True" , color = "r" , alpha = .5)[0]
np.histogram(data,bins,normed="True")
Out[269]:
In [ ]:
In [264]:
list_r = []
model = Zheng07(threshold = -21)
for i in xrange(100):
model.populate_mock()
gid = model.mock.compute_fof_group_ids()
x = np.histogram(richness(gid),bins , normed="True")[0]
list_r.append(x)
In [270]:
d_r = np.array(list_r)
print d_r.shape
plt.pcolor(d_r)
plt.colorbar()
Out[270]:
In [ ]:
In [591]:
plt.imshow(np.cov(d_r.T), interpolation = "None" , origin = "lower")
plt.colorbar()
Out[591]:
In [273]:
C_diag = np.diag(np.cov(d_r.T))+1.e-11
In [292]:
np.sum(((np.mean(d_r,axis=0)-y)**2./C_diag)[10:])
Out[292]:
In [291]:
np.sum(((np.mean(d_r,axis=0)-z)**2./C_diag)[10:])
Out[291]:
In [276]:
data_richness = np.mean(d_r,axis=0)
In [277]:
variance = C_diag
variance.shape , data_richness.shape
Out[277]:
In [278]:
covariance = np.cov(d_r.T)
In [279]:
np.savetxt("data_richness.dat" , data_richness)
In [280]:
np.savetxt("covariance_richness.dat" , covariance)
In [281]:
np.savetxt("variance_richness.dat" , variance)
In [ ]:
In [ ]:
In [293]:
In [ ]:
In [ ]:
In [301]:
zhengerror = np.array([3.46285183e+03, 1.64749125e+03, 7.80527281e+02,
3.30492728e+02, 1.38927882e+02, 5.91026616e+01,
2.45091664e+01, 1.10830722e+01, 5.76577829e+00,
3.14415063e+00, 1.88664838e+00, 1.07786531e+00,
5.54622962e-01, 2.87849970e-01])
In [508]:
covariance = np.cov(xr.T)
pylab.pcolor(covariance)
pylab.colorbar()
pylab.show()
icov = np.linalg.inv(covariance)
icov2 = np.linalg.solve(covariance , np.eye(14))
pylab.savetxt("inverse_clustering_covariance.dat" , icov2)
pylab.pcolor(icov)
pylab.colorbar()
pylab.show()
err = (np.diag(icov2))**-.5
print np.diag(covariance)
print err
In [316]:
m = Zheng07(threshold = -21)
In [319]:
m.param_dict["alpha"] = 2.
In [320]:
m.populate_mock()
In [321]:
r , x = m.mock.compute_galaxy_clustering()
rr , xx = model.mock.compute_galaxy_clustering()
In [609]:
print rr
pylab.figure(figsize(10,10))
pylab.errorbar(rr , xx , yerr = 2. * err , fmt="-ok", ms=1,
capsize=2, alpha=1.)
pylab.yscale("log")
pylab.xscale("log")
pylab.xlabel("$r$")
pylab.ylabel("$\\xi (r)$")
#pylab.xlim((,16))
Out[609]:
In [315]:
model.param_dict
Out[315]:
In [362]:
m = np.logspace(11,16, 100)
from scipy.special import erf
In [363]:
def mean_central(M):
logMmin = model.param_dict["logMmin"]
sigma = model.param_dict["sigma_logM"]
return .5 * (1. + erf((np.log10(M) - logMmin)/(sigma)))
In [ ]:
In [364]:
def mean_satelite(M):
ncen = mean_central(M)
M0 = 10. ** model.param_dict["logM0"]
M1 = 10. ** model.param_dict["logM1"]
alpha = model.param_dict["alpha"]
return ncen * ((M - M0)/(M1)) ** alpha
In [713]:
pylab.loglog(m , mean_central(m) , "r" , alpha = .1 , lw = 4.)
pylab.loglog(m , mean_satelite(m) , "b" , alpha = .1 , lw = 4.)
pylab.loglog(m , mean_central(m) + mean_satelite(m) , "g" , alpha = 1.)
pylab.ylabel("$N(M)$")
pylab.ylim((10**-6 , 10**10))
Out[713]:
In [351]:
print model.param_dict["logM0"]
In [451]:
model = Zheng07(threshold = -21)
In [453]:
model.populate_mock()
In [556]:
nbar = []
for i in xrange(500):
model.populate_mock()
nbar.append(model.mock.number_density)
In [557]:
nbar = np.array(nbar)
np.savetxt("mock_nbar.dat" , nbar)
In [558]:
print np.mean(nbar)
In [560]:
print np.var(nbar)
In [464]:
print model.param_dict
print model2.param_dict
In [ ]:
In [ ]:
In [515]:
np.diag(covariance)**.5
Out[515]:
In [516]:
np.savetxt("clustering_covariance.dat" , covariance)
In [518]:
model.mock.number_density
Out[518]:
In [526]:
np.logspace(np.log10(1e12),np.log10(10) , 30)
Out[526]:
In [ ]:
In [551]:
dist= [[ 9.34611722e+30 , 1.19139184e+03],
[ 4.40717814e+32 , 3.16556027e+03],
[ 6.97969451e+31 , 3.61021871e+02],
[ 6.23476568e+28 , 1.97302236e+04],
[ 3.05209675e+31 , 1.10247217e+03]]
eps= [ 1.00000000e+41 , 1.00000000e+12]
dist = np.array(dist)
In [553]:
[(dist[i]<eps).all() for i in range(dist.shape[0])]
Out[553]:
In [592]:
mean = np.mean(nbar)
var = np.var(nbar)
In [ ]:
In [568]:
In [570]:
In [571]:
mocks_xir = np.loadtxt("xir.dat")
data_xir = np.mean(mocks_xir , axis = 0)
covariance = np.loadtxt("clustering_covariance.dat")
cii = np.diag(covariance)
In [576]:
for i in range(200):
model.populate_mock()
n = model.mock.number_density
chi1.append((n - mean)**2. / var)
r , x = model.mock.compute_galaxy_clustering()
chi2.append(np.sum((data_xir - x)**2. / cii))
In [585]:
np.median(np.array(chi1))
Out[585]:
In [584]:
np.median(np.array(chi2))
Out[584]:
In [582]:
plt.hist(np.array(chi2))
Out[582]:
In [593]:
plt.hist(np.array(chi1))
Out[593]:
In [594]:
print "the mean of clustering distances with the correct HOD = " , np.array(chi2).mean()
print "the mean of number density distances with the correct HOD = " , np.array(chi1).mean()
In [ ]:
In [ ]:
In [605]:
model.populate_mock()
r , x = model.mock.compute_galaxy_clustering()
res = data_xir - x
%timeit np.dot(res , np.linalg.solve(covariance , res))
%timeit np.sum(res**2 / cii)
In [9]:
from halotools.mock_observables import wp
In [10]:
from halotools.mock_observables.clustering import tpcf
In [11]:
from halotools import mock_observables
In [ ]:
In [12]:
from halotools.empirical_models.mock_helpers import (three_dim_pos_bundle,
infer_mask_from_kwargs)
from halotools.mock_observables.clustering import wp
In [13]:
rbins = model_defaults.default_rbins
In [14]:
rbin_centers = (rbins[1:] + rbins[0:-1])/2.
In [15]:
print rbins
In [16]:
from halotools.sim_manager import supported_sims
In [17]:
cat = supported_sims.HaloCatalog()
In [18]:
L = cat.Lbox
mask = infer_mask_from_kwargs(model.mock.galaxy_table)
In [638]:
pos = three_dim_pos_bundle(table=model.mock.galaxy_table,
key1='x', key2='y', key3='z', mask=mask,
return_complement=False)
In [690]:
clustering = wp(pos, rbins, p_bins , period = L)
In [687]:
print clustering.shape
print rbin_centers.shape
print L
In [691]:
plot(rbin_centers , clustering)
plt.xscale("Log")
plt.yscale("Log")
In [689]:
p_bins = np.linspace(0,L/2,10)
In [694]:
pi_max = np.linspace(2, L/2 , 10)
cl = []
for pi in pi_max:
p_bins = np.linspace(0,pi,100)
clustering = wp(pos, rbins, p_bins , period = L)
cl.append(clustering)
In [750]:
figure = plt.figure(figsize = (10,10))
for i in range(0,10):
plot(rbin_centers , cl[i] , label = "$\pi_{max}$="+str(round(pi_max[i])))
#plot(rbin_centers , cl[i] , label = "$\pi_{max}$={0:.2f}".format(pi_max[i]))
plt.ylabel("$w_p$")
plt.xlabel("$r_p$")
plt.xscale("Log")
plt.yscale("Log")
plt.legend()
In [732]:
binsize = np.logspace(1,3,20)
cl2 = []
for bsize in binsize:
p_bins = np.linspace(0,125,bsize)
clustering = wp(pos, rbins, p_bins , period = L)
cl2.append(clustering)
In [754]:
figure = plt.figure(figsize = (10,10))
for i in range(0,10):
plot(rbin_centers , cl2[i] , label = "Bin size="+str(round(125/binsize[i]))+" Mpc")
#plot(rbin_centers , cl[i] , label = "$\pi_{max}$={0:.2f}".format(pi_max[i]))
plt.title("$\pi_{max}=$"+str(125)+" Mpc")
plt.ylabel("$w_p$")
plt.xlabel("$r_p$")
plt.xscale("Log")
plt.yscale("Log")
plt.legend()
In [739]:
binsize = np.logspace(1,3,20)
cl3 = []
for bsize in binsize:
p_bins = np.linspace(0,50,bsize)
clustering = wp(pos, rbins, p_bins , period = L)
cl3.append(clustering)
In [772]:
figure = plt.figure(figsize = (10,10))
for i in range(0,10):
plot(rbin_centers , cl3[i] ,label = "Bin size="+str(round(50/binsize[i]))+" Mpc")
#plot(rbin_centers , cl[i] , label = "$\pi_{max}$={0:.2f}".format(pi_max[i]))
plt.title("$\pi_{max}=$"+str(50)+" Mpc")
plt.ylabel("$w_p$")
plt.xlabel("$r_p$")
plt.xscale("Log")
plt.yscale("Log")
plt.legend()
In [786]:
model.populate_mock()
pos = three_dim_pos_bundle(table=model.mock.galaxy_table,
key1='x', key2='y', key3='z', mask=mask,
return_complement=False)
pbins = np.linspace(0 , 125 , 50)
clustering = wp(pos, rbins, p_bins , period = L)
In [770]:
from halotools.mock_observables import s_mu_tpcf
In [769]:
#s_bins = np.logspace(-2,-1,10)
#mu_bins = np.linspace(0,1,50)
#xis = s_mu_tpcf(pos, s_bins, mu_bins, period=L)
In [927]:
wp_list = []
p_bins = np.linspace(0,125,100)
for i in range(500):
model.populate_mock()
pos = three_dim_pos_bundle(table=model.mock.galaxy_table,
key1='x', key2='y', key3='z')
wp_list.append(wp(pos, rbins, p_bins , period = np.array([L,L,L])))
In [926]:
import time
a = time.time()
#pos = three_dim_pos_bundle(table=model.mock.galaxy_table,key1='x', key2='y', key3='z' ,return_complement=False)
#wp(pos, rbins, pi_bins , period = np.array([L,L,L]))
model.mock.compute_galaxy_clustering()
print time.time() - a
In [921]:
pi_bins = np.linspace(0,125,200)
d1 = wp(pos, rbins, pi_bins , period = np.array([L,L,L]))
In [922]:
pi_bins = np.linspace(0,125,100)
%timeit d2 = wp(pos, rbins, pi_bins , period = np.array([L,L,L]))
In [929]:
wps = np.array(wp_list)
In [931]:
wp_cov = np.cov(wps.T)
In [937]:
imshow(wp_cov , interpolation = "None" , origin = "lower")
colorbar()
Out[937]:
In [938]:
np.savetxt("wps.dat" , wps)
np.savetxt("wp_covariance.dat" , wp_cov)
In [946]:
data_wp = np.mean(wps , axis = 0)
cii = np.diag(wp_cov)
chi3 = []
for i in range(500):
#model.populate_mock()
#pos = three_dim_pos_bundle(table=model.mock.galaxy_table,
# key1='x', key2='y', key3='z')
#wp_list.append(wp(pos, rbins, p_bins , period = np.array([L,L,L])))
#r , x = model.mock.compute_galaxy_clustering()
chi3.append(np.sum((data_wp - wps[i])**2. / cii))
In [947]:
np.array(chi3).mean()
Out[947]:
In [957]:
model = Zheng07(threshold = -21)
model.populate_mock()
In [958]:
g = model.mock.compute_fof_group_ids()
In [959]:
r = richness(g)
In [979]:
print r.min()
In [980]:
bins = np.arange(1,46)
In [1017]:
histogram = plt.hist(r , bins)
plt.close()
number = histogram[0]
print histogram[1]
print number.shape
gmp = number/250**3.
figure = plt.figure(figsize = (10,10))
plt.semilogx(bins[:-1] , np.log10(gmp) , "bo-")
plt.xlabel("N")
plt.ylabel(r"$\log\, n_{grp}(N) \, [h^{3}Mpc^{-3}]$")
plt.show()
plt.close()
figure = plt.figure(figsize = (10,10))
plt.plot(bins[:-1] , np.log10(number) , "bo-")
plt.xlabel("N")
plt.ylabel(r"$\log\, Abundance(N)$")
plt.xlim((0,30))
plt.show()
plt.close()
In [1012]:
n_mocks = 20
grp_list = []
In [1305]:
print model.param_dict
#for k in range(n_mocks):
# model.populate_mock()
# gids = model.mock.compute_fof_group_ids()
# grp_list.append(richness(gids))
%timeit model.mock.compute_fof_group_ids(Nthreads = 4)
In [1070]:
#fig1, ax1 = plt.subplots()
print model.mock.number_density
In [1318]:
bins = [10,14,20,30]#[10,12,14,16,18,20,24,28]
gmps = []
for i in range(20):
r = grp_list[i]
histogram = plt.hist(r , bins)
#plt.show()
plt.close()
number = histogram[0]
print number
#print histogram[1]
#print number.shape
gmp = number/250**3.
#print np.sum(gmp) , gmp
gmps.append(gmp)
"""
#figure = plt.figure(figsize = (10,10))
fig1, ax1 = plt.subplots()
ax1.semilogx(bins[:-1] , number , "o-")
xticks = np.arange(10, 50, 10)
plt.xticks(xticks)
ax1.set_xscale('log')
ax1.set_xticks(xticks)
ax1.get_xaxis().set_major_formatter(matplotlib.ticker.ScalarFormatter())
ax1.set_xlim((10,50))
ax1.set_ylim((0,20))
ax1.set_xlabel("N")
ax1.set_ylabel("Abundance(N)")
ax1.set_title("Group multiplicity of mock "+str(i))
plt.show()
plt.close()"""
gmps = np.array(gmps)
print np.mean(gmps , axis = 0) , np.std(gmps , axis = 0)
In [ ]:
In [ ]:
In [1097]:
In [1099]:
In [1137]:
In [1108]:
In [1319]:
gmp_mean = gmps[0]
gmp_std = np.std(gmps , axis = 0)
In [1320]:
def distance(x):
return ((x - gmps[0])**2.) / gmp_std**2.
In [1321]:
for i in range(4):
model2 = Zheng07(threshold = -21)
model2.param_dict["alpha"] = 1.15
model2.populate_mock()
id2 = model2.mock.compute_fof_group_ids()
r2 = richness(id2)
n2 = plt.hist(r2 , bins)[0]/250.**3
plt.close()
print distance(n2) , np.sum(distance(n2))
In [1322]:
for i in range(4):
model2 = Zheng07(threshold = -21)
model2.param_dict["alpha"] = 1.3
model2.populate_mock()
id2 = model2.mock.compute_fof_group_ids()
r2 = richness(id2)
n2 = plt.hist(r2 , bins)[0]/250.**3
plt.close()
print distance(n2) , np.sum(distance(n2))
In [1323]:
for i in range(4):
model2 = Zheng07(threshold = -21)
model2.param_dict["alpha"] = 0.9
model2.populate_mock()
id2 = model2.mock.compute_fof_group_ids()
r2 = richness(id2)
n2 = plt.hist(r2 , bins)[0]/250**3
plt.close()
print distance(n2) , np.sum(distance(n2))
In [1293]:
for i in range(4):
model2 = Zheng07(threshold = -21)
model2.param_dict["alpha"] = 1.25
model2.populate_mock()
id2 = model2.mock.compute_fof_group_ids()
r2 = richness(id2)
n2 = plt.hist(r2 , bins)[0]/250**3
plt.close()
print distance(n2)
In [1324]:
m = Zheng07()
In [1325]:
m.param_dict
Out[1325]:
In [1326]:
m.populate_mock()
In [1327]:
m.threshold
Out[1327]:
In [1338]:
%timeit m.mock.compute_fof_group_ids()
In [1332]:
t = richness(g)
In [1340]:
plt.hist(t , np.arange(10,60 , 4))
Out[1340]:
In [1347]:
rich20 = []
for i in range(50):
m = Zheng07()
m.populate_mock()
gid = m.mock.compute_fof_group_ids()
rich = richness(gid)
rich20.append(rich)
#plt.hist(rich , np.arange(10,100))
#plt.xlim((10,100))
#plt.show()
In [1348]:
bins = np.arange(10,60,1)
gmp = []
for i in range(50):
gmp.append(plt.hist(rich20[i] , bins)[0]/250.**3.)
In [1351]:
plt.imshow(gmp , interpolation = "None")
plt.colorbar()
Out[1351]:
In [1357]:
for i in range(50):
plt.semilogx(bins[:-1] , np.log10(gmp[i]))
In [1440]:
sdss_bins = np.array([3,4,5,6,7,8,9,10,11,12,13,14,15,16,17,18,19,20,22,25,29,31,35,43,61])
sdss_bins = np.array([10,11,12,13,14,15,16,17,18,19,20,22,25,29,31,35,43,61])
In [1441]:
sdss_gmp = []
sdss_counts = []
for i in range(50):
sdss_gmp.append(plt.hist(rich20[i] , sdss_bins)[0]/250.**3.)
sdss_counts.append(plt.hist(rich20[i] , sdss_bins)[0])
plt.close()
In [1442]:
sdss_mean = np.mean(sdss_gmp , axis = 0)
sdss_count_mean = np.mean(sdss_counts ,axis = 0)
sdss_poisson = np.sqrt(sdss_count_mean) / 250. ** 3.
In [1443]:
sdss_sigma = np.std(sdss_gmp , axis =0)
print sdss_sigma , sdss_poisson
In [1445]:
plt.figure(figsize(10,10))
plt.errorbar(0.5*(sdss_bins[:-1]+sdss_bins[1:]) , sdss_mean, yerr = np.sqrt(sdss_poisson**2.+sdss_sigma**2.)
,fmt="ok", capsize=2.0)
plt.plot(0.5*(sdss_bins[:-1]+sdss_bins[1:]) , sdss_gmp[10] , alpha = .5)
plt.plot(0.5*(sdss_bins[:-1]+sdss_bins[1:]) , sdss_gmp[40] , alpha = .5)
plt.plot(0.5*(sdss_bins[:-1]+sdss_bins[1:]) , sdss_gmp[30] , alpha = .5)
plt.plot(0.5*(sdss_bins[:-1]+sdss_bins[1:]) , sdss_gmp[20] , alpha = .5)
plt.ylabel(r"Group Multiplicity Function [$h^{3}Mpc^{-3}$]")
plt.xlabel("Group Richness")
plt.xscale("log")
plt.yscale("log")
plt.xlim((10,60))
plt.title("Multiplicity function of the mock data with the true HOD")
plt.savefig("/home/mj/public_html/Multiplicity_function")
In [1449]:
def distance(model):
model.populate_mock()
group_ids = model.mock.compute_fof_group_ids()
R = richness(group_ids)
binned_R = plt.hist(R , sdss_bins)[0] / 250**3.
return np.sum((binned_R - sdss_mean)**2. / (sdss_sigma**2. + sdss_poisson**2.))
In [1463]:
q = Zheng07()
q.param_dict["alpha"] = 1.0
print distance(q)
In [1465]:
sdss_stdev = (sdss_sigma**2. + sdss_poisson**2.)**.5
sdss_bins = sdss_bins
sdss_mean = sdss_mean
np.savetxt("gmf_Mr20.dat" , sdss_mean)
np.savetxt("gmf_bins_Mr20.dat" , sdss_bins)
np.savetxt("gmf_noise_Mr20.dat" , sdss_stdev)
In [8]:
m = Zheng07()
print m.param_dict
In [ ]:
def ptpcf(tpcf , pbins , rbins):
In [ ]: