This notebook will make abundance matched catalogs for Jeremy and Zhongxu. I'm gonna send this notebook along to them as well in case there's something not quite right that they want to adjust. The catalogs requested were defined as follows:

2 catalogs, M_peak and V_max @ M_peak (which is how, I believe, V_max is defined in this catalog. Will check tho).
scatter 0.18 dex
number density of 4.2e-4
z = 0.55
using the SMF Jeremy provided, which is in this directory with name DR10_cBOSS_WISE_SMF_z0.45_0.60_M7.dat
On DS14, which is located on ki-ls at /nfs/slac/des/fs1/g/sims/yymao/ds14_b_sub, courtesy of Yao
Include in the catalog, along with the galaxies, M_vir, x, y, z, vx, vy, vz, M_gal, am_i_a_satellite?, and M_host



In [9]:

    
from os import path
import numpy as np
from AbundanceMatching import *
from halotools.sim_manager import RockstarHlistReader, CachedHaloCatalog



In [10]:

    
from matplotlib import pyplot as plt
%matplotlib inline
import seaborn as sns
sns.set()



In [11]:

    
#halo_dir = '/nfs/slac/des/fs1/g/sims/yymao/ds14_b_sub/hlists/'
halo_dir = '/scratch/users/swmclau2/hlists/ds_14_b_sub/hlists/'
a = 0.65
z = 1.0/a - 1 # ~ 0.55
fname = path.join(halo_dir,  'hlist_%.5f.list'%a)
print fname









    



/scratch/users/swmclau2/hlists/ds_14_b_sub/hlists/hlist_0.65000.list



In [12]:

    
columns_to_keep = {'halo_id': (1, 'i8'), 'halo_upid':(6,'i8'), 'halo_mvir':(10, 'f4'), 'halo_x':(17, 'f4'),\
                        'halo_y':(18,'f4'), 'halo_z':(19,'f4'),'halo_vx':(20,'f4'), 'halo_vy':(21, 'f4'), 'halo_vz':(22,'f4'),
                  'halo_rvir': (11, 'f4'),'halo_rs':(12,'f4'), 'halo_mpeak':(58, 'f4'),'halo_vmax@mpeak':(72, 'f4')}



In [13]:

    
columns_to_keep









    Out[13]:





{'halo_id': (1, 'i8'),
 'halo_mpeak': (58, 'f4'),
 'halo_mvir': (10, 'f4'),
 'halo_rs': (12, 'f4'),
 'halo_rvir': (11, 'f4'),
 'halo_upid': (6, 'i8'),
 'halo_vmax@mpeak': (72, 'f4'),
 'halo_vx': (20, 'f4'),
 'halo_vy': (21, 'f4'),
 'halo_vz': (22, 'f4'),
 'halo_x': (17, 'f4'),
 'halo_y': (18, 'f4'),
 'halo_z': (19, 'f4')}



In [14]:

    
simname = 'ds_14_b_sub'

Only run the below if you want to cache, which is useful maybe the first time (maybe). It takes ~30 min and some disk space, so be warned.

Update (Feb 1st, 2019): I had to edit halotools to make this work. The last line of the halocat was missing values... Specifically making the reader stop iteration once it encountered an indexerror.

reader = RockstarHlistReader(fname, columns_to_keep, '/scratch/users/swmclau2/halocats/hlist_%.2f.list.%s.hdf5'%(a, simname),\ simname,'rockstar', z, 'default', 1000.0, 2.44e9, overwrite=True, header_char = '#') reader.read_halocat(['halo_rvir', 'halo_rs'], write_to_disk=False, update_cache_log=False) reader.add_supplementary_halocat_columns() reader.write_to_disk() reader.update_cache_log()



In [15]:

    
halocat = CachedHaloCatalog(simname = simname, halo_finder='rockstar', redshift = z,version_name='default')

halocat.halo_table.colnames

%%bash pwd

np.min(halocat.halo_table['halo_mvir'])



In [18]:

    
smf = np.genfromtxt('DR10_cBOSS_WISE_SMF_z0.45_0.60_M7.dat', skip_header=True)[:,0:2]



In [19]:

    
smf[:5]









    Out[19]:





array([[  1.00000000e+01,   7.02900000e-08],
       [  1.01000000e+01,   3.11800000e-05],
       [  1.02000000e+01,   4.13400000e-06],
       [  1.03000000e+01,   1.22000000e-07],
       [  1.04000000e+01,   2.03000000e-05]])



In [20]:

    
nd = 4.2e-4 #nd of final cat



In [21]:

    
ab_property = 'halo_vmax@mpeak'



In [22]:

    
af = AbundanceFunction(smf[:,0], smf[:,1], (10, 12.9))

scatter = 0.18
remainder = af.deconvolute(scatter, 20)









    



/home/users/swmclau2/.local/lib/python2.7/site-packages/AbundanceMatching/AbundanceFunction.py:22: RuntimeWarning: overflow encountered in exp
  return -np.exp(a*x+b) + c*x + d



In [ ]:

    
nd_halos = calc_number_densities(halocat.halo_table[ab_property], 1000.0) #don't think this matters which one i choose here



In [ ]:

    
#check the abundance function
plt.semilogy(smf[:,0], smf[:,1], lw = 5)
x = np.linspace(10, 12.9, 101)
plt.semilogy(x, af(x))
plt.show()



In [ ]:

    
catalog = af.match(nd_halos, scatter)



In [ ]:

    
catalog.shape



In [ ]:

    
n_obj_needed = int(nd*(1000.0**3))



In [ ]:

    
sort_idxs = np.argsort(catalog)
final_catalog = catalog#[sort_idxs[:n_obj_needed]]

plt.hist(final_catalog, bins = 500);



In [ ]:

    
output = halocat.halo_table#[sort_idxs[:n_obj_needed]]



In [ ]:

    
output['gal_smass'] = final_catalog



In [ ]:

    
np.sum(np.isnan(final_catalog))/(1.0*len(final_catalog))

output.write('catalog_ab_%s.hdf5'%ab_property, format = 'hdf5', path = './catalog_ab_%s.hdf5'%ab_property, overwrite=True)



In [ ]:

    
print ab_property



In [ ]: