-
Notifications
You must be signed in to change notification settings - Fork 0
/
082117_wp_ggl_uniform_random_AB.py
172 lines (116 loc) · 5.15 KB
/
082117_wp_ggl_uniform_random_AB.py
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
144
145
146
147
148
149
150
151
152
153
154
155
156
157
158
159
160
161
162
163
164
165
166
167
168
169
170
171
172
import argparse
parser = argparse.ArgumentParser(description='Compute delta sigma')
parser.add_argument('--td',type=int,default=0,dest='td',help='time delay in seconds')
parser.add_argument('--lower',required=True,nargs=7,type=float,dest='lower')
parser.add_argument('--upper',required=True,nargs=7,type=float,dest='upper')
parser.add_argument('--outfile',required=True,dest='outfile')
args = parser.parse_args()
import time
time.sleep(args.td)
import sys
import collections
import gc
import numpy as np
from concurrent.futures import ProcessPoolExecutor as Pool
from halotools.sim_manager import CachedHaloCatalog
from HOD_models import decorated_hod_model
from HOD_models import standard_hod_model
from halotools.empirical_models import MockFactory
from halotools.mock_observables import return_xyz_formatted_array
from halotools.mock_observables import delta_sigma
from halotools.mock_observables import wp
from halotools.utils import randomly_downsample_data
from halotools.utils import crossmatch
from Ngal_estimate import Ngal_estimate
##########################################################
param_names = ('alpha','logM1','sigma_logM','logM0','logMmin','mean_occupation_centrals_assembias_param1','mean_occupation_satellites_assembias_param1')
output_names = ('ngals','deltasigma','rp','wprp','param')
##########################################################
halocat = CachedHaloCatalog(simname = 'diemerL0500', version_name = 'antonio', redshift = 0, \
halo_finder = 'rockstar',ptcl_version_name='antonioz0')
Lbox = 500
particle_portion = 0.1
rp_bins = np.logspace(-1.398, 1.176, 14) ##to match the leauthaud paper
num_ptcls_to_use = int(1e6)
particle_masses = np.zeros(num_ptcls_to_use)+halocat.particle_mass/particle_portion
total_num_ptcls_in_snapshot = len(halocat.ptcl_table)
downsampling_factor = total_num_ptcls_in_snapshot/float(num_ptcls_to_use)
##ggl
pi_max = 60
r_wp = np.logspace(-1, np.log10(Lbox)-1, 20)
##wp
pos_part = return_xyz_formatted_array(*(halocat.ptcl_table[ax] for ax in 'xyz'), period=Lbox)
pos_part = randomly_downsample_data(pos_part, num_ptcls_to_use)
#########################################################
def calc_all_observables(param):
model.param_dict.update(dict(zip(param_names, param))) ##update model.param_dict with pairs (param_names:params)
c = Ngal_estimate(halocat,param)
n_est = c.ngal_estimate()
if n_est<1e5 and n_est>7.8e4:
try:
model.mock.populate()
except:
model.populate_mock(halocat)
gc.collect()
output = []
if model.mock.galaxy_table['x'].size<9.8e4 and model.mock.galaxy_table['x'].size>8e4:
pos_gals = return_xyz_formatted_array(*(model.mock.galaxy_table[ax] for ax in 'xyz'),period=Lbox)
pos_gals = np.array(pos_gals,dtype=float)
pos_gals_d = return_xyz_formatted_array(*(model.mock.galaxy_table[ax] for ax in 'xyz'), \
velocity=model.mock.galaxy_table['vz'], velocity_distortion_dimension='z',\
period=Lbox) ##redshift space distorted
pos_gals_d = np.array(pos_gals_d,dtype=float)
#ngals
output.append(model.mock.galaxy_table['x'].size)
#delta sigma
deltasigma = delta_sigma(pos_gals, pos_part, particle_masses=particle_masses, downsampling_factor=downsampling_factor, rp_bins=rp_bins, period=Lbox)
output.append(deltasigma[1])
output.append(deltasigma[0])
# wprp
output.append(wp(pos_gals_d, r_wp, pi_max, period=Lbox))
# parameter set
output.append(param)
else:
output = 1
return output
############################################################
def main(output_fname):
"""
default_priors['alpha']=[0.0,2.0]
default_priors['logM1']=[10.7,15.0]
default_priors['sigma_logM']=[0.02,1.5]
default_priors['logM0']=[9.0,14.0]
default_priors['logMmin']=[9.0,14.0]
default_priors['mean_occupation_centrals_assembias_param1']=[-1.0,1.0]
default_priors['mean_occupation_satellites_assembias_param1']=[-1.0,1.0]
"""
lower = np.array(args.lower[:])
upper = np.array(args.upper[:])
# p_step = (upper-lower)/float(args.nstep-1)
# grid = np.mgrid[0:args.nstep,0:args.nstep,0:args.nstep,0:args.nstep,0:args.nstep,0:args.nstep,0:args.nstep]
# grid = grid.reshape(7,args.nstep**7)
# grid = grid.T
# params = lower+grid*p_step
p_range = upper-lower
dev = np.random.rand(10000,7)
params = lower+dev*p_range
nproc = 50
global model
model = decorated_hod_model()
output_dict = collections.defaultdict(list)
with Pool(nproc) as pool:
for i, output_data in enumerate(pool.map(calc_all_observables, params)):
if i%50 == 0:
sys.stdout.write("\r%i" % i)
sys.stdout.flush()
if output_data!=1:
for name, data in zip(output_names, output_data):
output_dict[name].append(data)
for name in output_names:
output_dict[name] = np.array(output_dict[name])
np.savez(output_fname, **output_dict)
if __name__ == '__main__':
main(args.outfile)
with open(args.outfile+'_log','w') as f:
for arg in vars(args):
f.write(str(arg)+':'+str(getattr(args, arg))+'\n')